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I.  Introduction 


The  ultimate  goal  of  this  project  was  to  develop  a  device  for  encoding  and  sending 
modulated  high  frequency  optical  signals  utilizing  a  novel  combination  of  chemical  oscillators, 
chemiluminescent  reactions  and  microfluidic  technology.  Our  research  team  consisted  of 
chemists,  physicists  and  computer  scientists  from  Brandeis  University,  supported  by  scientists 
from  Spectral  Sciences,  Inc.  (SSI)  and  the  Air  Force  Research  Laboratory,  Hanscom  AFB,  as 
well  as  engineers  from  RainDance  Technologies,  Inc.  (RDT). 

The  problem  is  an  extremely  challenging  one  and  requires  the  solution  of  several 
subproblems  along  with  integration  of  those  solutions  to  produce  a  viable  device.  Among  the 
tasks  that  need  to  be  accomplished  are:  identification  and  testing  of  an  appropriate  chemical 
reaction  or  set  of  reactions  capable  of  producing  sufficiently  intense  chemiluminescence  in  an 
appropriate  wavelength  range  and  with  a  useful  energy  capacity  (at  least  100  gW/g);  design  and 
testing  of  a  microfluidic  circuit  that  makes  it  possible  to  generate  the  desired  signal;  development 
of  a  system  that  permits  the  encoding  and  repeated  transmission  of  messages  of  60  characters  or 
more  from  an  alphanumeric  character  set;  construction  (or  identification)  and  integration  into 
the  circuit  of  a  chemical  power  source  to  propel  the  luminescent  material  around  the  circuit; 
development,  analysis  and  modeling  of  sensors  and  algorithms  optimized  for  detection  of  the 
signals  generated  by  the  device,  including  construction  of  a  programmable,  radiometrically 
calibrated,  modulated  light  source  that  is  capable  of  emulating  the  proposed  sources  and 
encoding  concepts. 

We  have  made  considerable  progress  toward  solving  several  of  these  problems  and  have 
discovered  several  potentially  significant  new  phenomena  in  the  course  of  our  investigations. 
We  report  these  developments  below.  Nonetheless,  there  remains  significant  work  to  be  done  in 
order  to  produce  an  operable  device.  For  each  aspect  of  the  problem,  we  outline  below  the 
approaches  taken  and  the  results  obtained,  along  with  brief  conclusions  about  prospects  for 
future  success. 

II.  Research  Accomplished 

A.  Chemiluminescent  Reactions 

A  key  element  of  our  approach  is  the  generation  of  light  by  a  chemiluminescent  chemical 
reaction.  We  identified  and  tested  several  such  reactions,  both  in  the  near  infrared,  the 
wavelength  region  that  we  ultimately  hope  to  exploit,  and  in  the  visible,  which  is  more 
convenient  for  testing  the  microfluidic  circuit  and  other  elements  of  our  device,  and  which  may 
ultimately  prove  useful  in  itself. 

1.  Singlet  oxygen  generation 

The  most  promising  route  to  generating  continuous  high  intensity  chemiluminescence  in 
the  near  infrared  region  (NIR)  of  the  spectrum  is  through  the  production  of  singlet  oxygen.  The 
radiative  decay  of  singlet  O2  ('Ag  — »3Zg')  results  in  emission  in  the  NIR  at  A  =  1270  nm,  a 
wavelength  that  is  ideal  for  detection.  In  CCI4,  the  quantum  yield  of  this  NIR  phosphorescence 
is  about  0.07,  though  this  solvent  is  not  ideal  because  of  its  toxicity.  The  thermal  decomposition 
of  endoperoxides  appears  to  be  the  most  convenient  approach  to  generating  singlet  oxygen  in 
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situ,  and  we  succeeded  in  synthesizing  300  mg  of  1 ,4-dimethylnaphthalene  endoperoxide 
(DMNE),  the  most  promising  compound  we  identified. 

2.  Glowstick  Technology 

A  second  approach  to  producing  a  continuous  source  of  chemiluminescence  is  the  technology 
utilized  in  the  production  of  commercial  glowsticks.  Visible  light  of  various  colors  can  be 
generated  by  the  reaction  of  hydrogen  peroxide  with  a  phenyl  oxalate  ester  in  the  presence  of  a 
suitable  fluorescent  dye.  The  scheme  is  shown  below. 


We  purchased  red  glowsticks  from  a  commercial  vendor  and  extracted  the  reagents  to 
give  a  convenient  and  inexpensive  source  of  chemiluminescent  material  for  testing  our 
microfluidic  configurations.  We  also  discovered  a  vendor  who  manufactures  infrared 
glowsticks  for  the  military  (http://glowproducts.eom/glowsticks/infraredglowsticks/T  We 
acquired  a  stock  of  these  and  ran  some  preliminary  tests,  which  suggested  that  the  quantum  yield 
is  somewhat  less  than  that  of  the  chemiluminescence  produced  by  DMNE. 

3.  Luminescent  oscillators 

An  intriguing  possibility  is  to  utilize  a  chemical  oscillator  as  the  source  of 
chemiluminescence.  We  have  developed  a  chemiluminescent  oscillator  by  coupling  luminol  to 
the  Cu(II)-catalyzed  reaction  of  hydrogen  peroxide  with  potassium  thiocyanate  in  sodium 
hydroxide  (S.  Sattar  and  I.  R.  Epstein,  J.  Phys.  Chem.  94,  275-277  (1990)).  The  major 
limitations  on  employing  this  reaction  are  that  it  requires  a  continuous  flow  of  reactants  to 
sustain  the  oscillation  and  it  oscillates  relatively  slowly,  with  a  period  of  many  minutes.  Recent 
work  in  Iran  has  utilized  a  variety  of  organic  bases  in  place  of  NaOH,  resulting  in  a  significant 
increase  in  frequency  (A.  Samadi-Maybodi  and  R.  Akhoondi, 
Luminescence  23  ,  42  (2007))  as  well  as  in  intensity  and  duration  of  oscillation  in  a  closed 
system  (M.H.  Sorouraddin,  M.  Iranifam  and  K.  Amini,  J.  Fluor.  ]_8,  443  (2008).  We  carried  out 
several  initial  experiments  to  investigate  the  feasibility  of  producing  a  fast,  long-lasting 
chemiluminescent  oscillator,  but  to  date  have  not  succeded  in  increasing  the  lifetime  enough  to 
make  this  approach  feasible. 
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4.  Hydrogels 

A  potentially  more  attractive  approach  involves  coupling  the  long-lasting  BZ  oscillator  to  the 
luminol  reaction  in  a  supramolecular  hydrogel.  The  idea  is  to  utilize  the  periodic  oxidation  and 
reduction  of  the  BZ  catalyst  to  periodically  switch  the  chemiluminescence  reaction  on  and  off  in 
order  to  create  a  specific  pattern  of  light  emission  in  the  hydrogel.  Remarkably,  the  exceptional 
stability  and  unique  nanofiber  matrix  of  the  hydrogel  not  only  give  it  stability  to  the  extreme 
conditions  (pH  0-1)  of  the  BZ  system,  but  also  enhance  the  intensity  of  chemiluminescence  by 
improving  quantum  yield,  extend  the  lifetime  of  both  the  catalysts  of  chemiluminescence  and  the 
oscillating  reactions  by  minimizing  side  reactions,  and  allow  us  to  control  the  pattern  of  emission 
by  tailoring  the  molecular  hydrogelators.  We  can  maximize  the  compatibility  of  the  chemical 
oscillating  systems  with  devices  by  tuning  the  rheological  properties  of  the  hydrogels. 

The  Scheme  at  the  right  outlines  the 
design  principle  for  generating  the  oscillating 
chemiluminescent  supramolecular  hydrogel,  a 
heterogeneous  system  consisting  of  a 
nanofiber  matrix  and  a  liquid  phase 
immobilized  by  the  matrix.  The  substrate  and 
the  catalyst  for  the  chemiluminescent  reaction 
can  be  localized  in  the  nanofiber,  while  the 
oscillating  reaction  occurs  in  the  liquid  phase. 

Because  of  the  nanosize  diameter  of  the  fibers  and  the  liquid  phase,  the  substrate,  the  catalyst, 
and  the  components  of  the  oscillating  reaction  essentially  behave  as  if  they  were  on  a  surface  but 
able  to  interact  with  each  other.  Such  organization  mimics  the  environment  of  bioluminescence, 
enabling  us  to  maximize  the  quantum  yield  of  chemiluminescence  and  to  minimize  the  decay  of 
the  oscillating  reaction  by  confining  the  BZ  reactants  to  the  liquid  phase.  To  tune  the  wavelength 
of  the  emission  into  the  desired  range,  we  can  use  the  chemiluminescence  to  excite  quantum  dots 
incorporated  into  the  nanofibers. 

5.  Conclusion 

If  IR  luminescence  is  essential,  the  DMNE  system  is  probably  the  best  option.  Other  less 
toxic  aprotic  solvents  should  be  investigated  to  avoid  the  problems  with  CCI4.  If  an  oscillatory 
signal  is  desirable,  it  might  be  possible  to  couple  this  reaction,  perhaps  in  a  gel  system,  to  a 
chemical  oscillator.  For  visible  emission,  luminol-based  reactions  offer  the  advantage  of  good 
quantum  yield  and  existing  oscillatory  behavior.  Use  of  a  hydrogel  appears  to  offer  the  most 
promising  path  to  increasing  the  lifetime  of  the  oscillations  and  the  quantum  yield  of  emission. 


B.  Fast  Chemical  Oscillation 

The  Belousov-Zhabotinsky  (BZ)  reaction,  which  involves  the  oxidation  by  bromate  of  an 
organic  substrate,  typically  malonic  acid,  in  the  presence  of  a  metal  ion  or  complex  as  catalyst  in 
a  concentrated  solution  of  sulfuric  acid,  is  the  best  known  and  most  versatile  and  robust  of  the 
oscillating  chemical  reactions.  Unlike  most  chemical  oscillators,  which  require  a  continuous 
flow  of  fresh  reactants  to  sustain  oscillatory  behavior  for  more  than  a  few  cycles,  the  BZ  system 
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will  oscillate  for  several  hours  and  hundreds  of  cycles  in  a  closed  vessel  under  appropriate 
conditions.  On  the  other  hand,  typical  BZ  oscillations  have  periods  of  several  minutes  or  longer, 
i.e.,  frequencies  of  several  millihertz.  Since  we  sought  to  generate  signals  with  frequencies  in  the 
range  of  10  Hz  or  more,  we  decided  to  try  to  accelerate  the  BZ  oscillations  by  three  to  four 
orders  of  magnitude.  The  key  elements  in  our  scheme  were  to  work  at  higher  than  normal  acid 
concentrations  and  temperatures  and  to  vary  the  concentrations  of  the  reactants  and  catalyst  as 
well  as  the  identity  of  the  catalyst. 

1.  Experiments 

We  carried  out  an  extensive  investigation  of  the  parameter  space  of  the  reaction,  varying 
temperature,  acidity,  identity  of  the  catalyst,  identity  of  the  substrate,  and  reactant 
concentrations.  We  explored  the  use  of  1,4-cyclohexanedione  (CHD)  and  acetylacetone  as 
alternatives  to  the  usual  malonic  acid  (MA)  substrate,  but  found  that  malonic  acid  gives  the  best 
performance.  Experiments  were  conducted  with  cerium  (III),  manganese  (II)  and  ferroin  as 
catalysts.  Temperature  was  varied  from  room  temperature  to  over  70  °C,  above  which 
oscillations  ceased  or  became  difficult  to  monitor  because  of  the  production  of  CO2  bubbles. 

Table  1  illustrates  typical  data  with  CHD  as  substrate  and  ferroin  as  catalyst.  Our  best  results 
were  obtained  with  MA  and  ferroin  at  a  temperature  of  60  °C,  where  a  frequency  of  1 8  Hz  was 
observed. 

Table  1.  Oscillations  in  the  ferroin-catalyzed  BZ-CHD  system 


Initial  Concentrations  (M) 

T(°C) 

F  (Hz) 

NaBr03 

CHD 

Ferroin 

h2so4 

0.5 

0.6 

0.001 

1 

40 

0.5 

1 

1.2 

0.001 

1 

40 

no  oscillation 

0.5 

0.6 

0.0005 

2 

40 

1.5 

0.4 

0.6 

0.0005 

2 

40 

1.2 

0.3 

0.6 

0.0005 

2 

40 

1 

0.4 

0.6 

0.0005 

3 

40 

2 

0.4 

0.6 

0.0005 

4 

40 

no  oscillation 

0.4 

0.6 

0.0005 

2 

50 

4.5 

0.4 

0.6 

0.0005 

2 

60 

6 

0.4 

0.6 

0.0005 

2 

65 

10 

0.4 

0.6 

0.0005 

2 

70 

13 

0.4 

0.625 

0.0005 

2 

70 

12 

0.425 

0.6 

0.0005 

2 

70 

12 

0.5 

0.6 

0.0005 

2 

70 

no  oscillation 

Figure  1  shows  an  example  of  high  frequency  oscillations.  We  were  indeed  able  to  increase 
the  frequency  of  a  BZ  oscillator  by  roughly  a  factor  of  5,000,  but  20  Hz  appears  to  represent  the 
practical  upper  limit  with  this  approach,  and  in  a  closed  system  these  high  frequency  oscillations 
cannot  be  sustained  for  more  than  a  few  minutes,  though  with  a  flow  they  can  presumably 
continue  indefinitely. 
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Figure  1.  High  frequency  (18  Hz)  oscillations  in  the  ferroin-catalyzed  BZ  system.  Intiital 
concentrations:  [MA]0  =  0.5  M;  [NaBrO3]0=  0.5  M;  [H2SC>4]o  =  1.14  M;  [Ferroin]  =  0.0025  M.  T  = 
60  °C.  Noisy  data  result  from  production  of  CO  bubbles,  which  can  be  reduced  by  adding  a  surfactant. 

2.  Simulations 

We  developed  a  model  incorporating  all  reactant  species  based  on  earlier  work  by 
Zhabotinsky  et  al.  to  describe  the  shape  of  oscillations  and  wave  properties  in  the  BZ  reaction  at 
elevated  temperatures.  The  simplified  set  of  equations  used  to  describe  our  experimental  findings 
is 


cU7clt=  -k2h  AY+k3h0A  Y-2luX2-k5h0AX+k.5  U2+  k6U(C-Z)-k.(yXZ  (1) 

di7dt=  -k2hoAT-k3hoA7+qk7BZ+k9B  (2) 

dZ/dt=k6  £/(C-Z)-k_6XZ-k7BZ  (3) 

d  U! dt=2k5hoAX-2k.5  U2-k6U(C-Z)+k.6XZ  (4) 


where  X=  [HBr02],  Y  =  [Br  ],  Z  =  [Ce4+],  U  =  [HBr02],  h0  is  Hammet’s  acidity  function,  A 
=  [NaBr03]  =  h0/(0.2+h0)  [NaBrO3]0,  B.  =  [MA],  C  =  [Ce3+]  +  [Ce4+],  and  k4*  =  k4(l+0.87h0). 
For  numerical  calculations,  based  on  the  initial  reactant  concentrations,  we  assigned  A0=  0.5  M, 
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Bo  =  2  M,  C  =  O.OOl  M,  ho  =  3,  which  corresponds  to  [H2S04]=1.5  M.  The  stoichiometric 
parameter  q  was  set  to  0.6. 

Activation  energies  calculated  by  fitting  our  data  were  in  good  agreement  with  data 
available  in  the  literature,  and  the  simulations  gave  good  agreement  with  experiments  on  the 
ferro in-catalyzed  malonic  acid  system,  as  illustrated  in  Figure  2. 


0.05  ^ 

© 

~o 

0.025  3 

Q. 

E 

o  * 


40  50  60  70  80  90 

temperature  (°C) 


Figure  2.  Experimental  (left)  and  simulated  (right)  frequency  (squares)  and  amplitude 
(triangles)  of  oscillation  in  the  BZ  system  as  a  function  of  temperature. 

3.  Conclusion 

We  have  succeeded  in  increasing  the  frequency  of  oscillation  above  10  Hz  by  manipulating 
the  concentrations  and  temperature.  Under  these  conditions,  the  practical  upper  limit  of 
frequency  is  about  20  Hz,  and  the  lifetime  of  the  oscillation  decreases  significantly.  The 
simulations  suggest  that  higher  frequencies  might  be  obtainable  by  putting  the  system  under 
pressure,  which  would  allow  the  temperature  to  be  increased  beyond  100  °C.  Extending  the 
lifetime  would  require  use  of  a  flow  system. 


C.  Microfluidic  Circuitry 

A  key  aspect  of  our  approach  was  to  utilize  microfluidic  techniques  to  produce  controlled 
flows  of  chemiluminescent  solutions  that  will  be  used  to  encode  the  signals.  In  collaboration 
with  engineers  at  RDT,  we  developed  and  tested  several  circuits  to  generate  and  propel 
continuous  streams  of  solution  or  streams  consisting  of  aqueous  droplets  separated  by  oil.  Since 
the  BZ  reactants  are  polar,  the  reaction  occurs  only  in  the  aqueous  droplets,  though 
communication  is  possible  between  droplets  via  nonpolar  intermediates  such  as  Br2  that  are 
generated  in  the  reaction  and  then  partition  into  the  oil  phase.  Since  the  droplet  diameter  is 
comparable  to  the  diffusion  length  (~  100  pm),  the  droplets  are  internally  synchronized  by 
diffusion  and  require  no  stirring. 


Final  Report  W911NF-07-1-0639 


Page  9 


1.  One-dimensional  experiments 

The  simplest  configuration  we  explored,  and  the  one  that  lends  itself  most  easily  to  the 
straightforward  approach  of  flowing  a  stream  of  drops  behind  a  coded  mask,  consists  of  a  narrow 
tube,  i.e.,  a  pseudo-one-dimensional  (ID)  arrangement,  which  is  filled  using  microfluidic 
techniques  with  a  fluid  composed  of  aqueous  drops  containing  the  luminescent,  oscillatory 
solution  separated  by  an  oil  phase.  The  fluid  moves  at  controllable  velocities  of  the  order  of 
cm/sec,  and  drop  sizes  can  easily  be  varied  in  the  50-300  pm  range.  The  circuits  are  constructed 
from  polydimethylsiloxane  (PDMS). 

Figure  3  shows  examples  of  some  tubes  containing  arrays  of  BZ  droplets.  Figure  4  shows 
the  droplets  being  generated.  Table  2  gives  a  sense  of  the  range  of  conditions  explored  and  url’s 
for  videos  of  some  of  the  experiments. 


Figure  3.  Tubes  containing  aqueous  droplets  of  varying  size  and  inter-droplet  distance  with  BZ  reagent 
(dark,  convex  objects)  separated  by  oil  (light,  concave).  Length  of  tube  shown  is  4.1  mm;  inner  diameter 
of  each  tube  is  150  pm. 


pp8 . 93o3 . 85wpsi_l  5f  ps- 1  knf 


Figure  4.  Generation  of  droplets  in  the  microfludic  circuit. 
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Table  2.  A  sampling  of  microfluidics  experiments 


Fluid  1 
for  drop 
genera¬ 
tion 

Fluid2 

Droplet 

Fluid3  for 
accelera¬ 
tion 

flow-rate 
(ul/hr),  pressure 
(psi) 

fluid  l/fluid2/ 
fluid3 

Frequency  (Hz) 
Gen/accel 
s=stable, 
a=aperiodic 

Speed 

(mm/s) 

Size 

(microns) 

video 

Oil 

Water 

Oil  (flow) 

250/100/500 

si  .25/1 .25 

1.75/4.5 

500x335 

(flow) 

(flow) 

250/100/1000 

si. 2/1. 2 

1.75/6.6 

500x335 

http://people.brandeis.edu/~ochoa/DARP 

A/250ol-lhw-lko2  30fbs-150nf.AVI 

250/100/2000 

si. 8/1. 8 

1.75/13 

500x335 

250/100/4000 

si. 66/1.66 

1.75/44 

500x335 

250/100/7000 

al. 76/1. 76 

1 .75/ — 5 
5 

500x335 

http://people.brandeis.edu/~ochoa/DARP 

A/250ol-lhw-7ko2  lkfps-3knf.AVI 

500/300/10000 

s2. 75/2. 75 

2.9/43 

500x335 

500/500/20000 

s3. 75/3. 75 

- 

545x335 

2000/2000/400 

00 

S14.25/14.25 

- 

548x335 

http://people.brandeis.edu/~ochoa/DARP 

A/2ko  1  -2kw-40ko2  lkfps-4knf.AVI 

5000/500/2000 

S12.2/12.2 

25.6/ 

293 

http://people.brandeis.edu/~ochoa/DARP 

A/5kol  5hw  2ko2  lkfps  lknf.AVI 

10000/1000/20 

00 

s25/25 

265 

15000/1500/20 

00 

s45/45 

240 

20000/2000/50 

00 

s60/60 

240 

30000/3000/50 

00 

s99/99 

223 

http://people.brandeis.edu/~ochoa/DARP 

A/3 0k  3k  5k  lkfps  lknf.AVI 

Oil 

(flow) 

Air  (flow) 

200/80/2000 

1.3/-1.3  unstable 
alternating  cycle 

520x335 

http://people.brandeis.edu/~ochoa/DARP 

A/2ho-8dw-2ka  30fbs-1800nf  l.AVI 

200/80/5000 

-/~1.4  very 
irregular 

Air 

(pressure) 

250/80/2psi 

ul.32/aperiodic 
drop  breaking 

530x335 

Oil 

Water 

- 

2.76/1.4 

a0.345 

470x335 

http://people.brandeis.edu/~ochoa/DARP 

A/pp2.76ol.4wpsi  15fps-lknf.AVI 

(pressur 

(pressur 

3 .2/2.3 

a0.286 

-470x335 

e) 

e) 

3.2/2.36 

s0.42 

-536x335 

3. 5/2.3 

a0.58 

355x335 

3.52/1.43 

a0.495 

357x335 

http://people.brandeis.edu/~ochoa/DARP 

A/pp3.52ol.43wpsi  15fps-lknf.AVI 

5/2.3 

al .  1 8 

305 

5.33/3.85 

a*0.588 

625x335 

6.25/2.97 

a*  1.8 

293 

6.36/3.86 

s3.1 

424x335 

6.86/5.64 

a*0.416 

530x335 

http://people.brandeis.edu/~ochoa/DARP 

A/pp6.86o5.64wpsi  10fps-6hnf.AVI 

7.93/5.64 

a*2.4 

514x335 
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Perhaps  surprisingly,  the  dominant  behavior  for  droplets  separated  by  less  than  400  pm  is 
out-of-phase  coupling.  While  intuitively  one  might  expect  that  bringing  the  droplets  close 
together  would  produce  in-phase  behavior,  this  turns  out  not  to  be  the  case  in  the  present 
configuration,  at  least  down  to  separations  of  50  pm.  Apparently,  the  fact  that  the  coupling 
occurs  via  bromine,  an  inhibitor  of  the  reaction,  is  responsible  for  this  unanticipated  behavior. 
Figure  5  shows  some  typical  results. 


ujji  ii»  ijj  piimiihb 


Figure  5.  Space-time  plots  showing  (a)  anti-phase  oscillations  with  spikes  of  oxidation  of  ferroin  seen  as 
light  horizontal  lines  across  BZ-droplets.  (b)  stationary  Turing  structure  with  alternating  oxidized  and 
reduced  states  evolving  from  initial  oscillatory  state.  Horizontal  length  of  the  frame  and  the  capillary  ID 
are  4.8  mm  and  150  pm,  respectively;  the  total  times  shown  in  (a)  and  (b)  are  5200  s  and  10800  s, 
respectively.  Patterns  extend  to  the  left  and  right  of  the  segments  shown. 

In  the  figure,  the  droplets  are  the  darker  vertical  columns,  separated  by  the  lighter  oil 
columns.  The  light  narrow  horizontal  stripes  are  the  brief  flashes  of  oxidized  ferroin  (ferriin) 
during  each  cycle.  These  spikes  start  off  in-phase  because  the  system  is  initiated  with  a  flash  of 
light  to  synchronize  the  drops.  In  frame  (b),  where  the  concentrations  are  slightly  different  from 
frame  (a),  the  oscillations  cease  after  some  time,  and  the  system  enters  a  stationary,  but 
patterned,  “Turing”  state,  where  alternate  drops  are  either  in  the  oxidized  or  the  reduced  state. 

2.  Two-dimensional  experiments 

We  also  used  microfluidic  techniques  to  prepare  two-dimensional  (2D)  arrays  of  BZ  droplets 
in  the  fluorinated  oil  HFE-7500.  Of  course,  2D  arrays  can  give  rise  to  a  much  richer  set  of 
phenomena  than  ID  arrays,  and  we  have  only  begun  to  scratch  the  surface  of  what  is  possible. 
Earlier  experiments  at  Brandeis  on  microemulsions  consisting  of  nanometer  diameter  droplets  of 
a  photosensitive  aqueous  BZ  solution  in  oil  showed  that  this  medium  can  serve  as  an  erasable 
memory  device  (A.  Kaminaga,  V.  K.  Vanag,  and  I.  R.  Epstein,  Angew.  Chem.  Int.  Ed.  45,  3087 
(2006)).  In  those  experiments,  however,  the  droplets  undergo  Brownian  motion.  Here,  with  the 
droplets  having  diameters  of  roughly  100  pm,  the  droplets  are  stationary,  and  it  should  be  easier 
to  control  the  patterns.  We  observed  coherent  behavior,  including  regions  of  in-  and  out-of- 
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phase  oscillation  and  traveling  waves.  Figure  6  shows  an  example.  In  other  experiments,  we 
observed,  after  the  droplets  spontaneously  self-assembled  into  a  hexagonal  array,  a  variety  of 
coherent  patterns,  the  most  interesting  of  which  consisted  of  a)  a  three-phase  oscillating  pattern, 
shown  schematically  in  Figure  7,  in  which  all  droplets  oscillate,  with  nearest  neighbors  being 
120°  out-of-phase  with  one  another,  and  b)  a  two-phase  mixed  pattern,  in  which  steady  state 
droplets  alternate  with  oscillatory  droplets  in  such  a  way  that  neighboring  oscillatory  droplets  are 
180°  out-of-phase,  and  the  pattern  is  symmetric  to  rotations  by  120°. 


Figure  6.  Traveling  wave  (white  area  at  lower  left)  of  oxidation  in  a  2D  array  of  BZ  water  droplets. 
White  “stars”  and  large  white  circle  are  defects  and  air  bubbles. 


Figure  7.  Schematic  drawing  of  the  three-phase  pattern. 
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In  addition  to  the  emulsion- like  arrangement  shown  in  Figure  6,  we  also  explored 
microfluidic  configurations  in  which  droplets  can  be  placed  and  then  connected  via  narrow 
channels  in  any  desired  2D  arrangement.  One  such  arrangement  is  shown  in  Figure  8. 


600  pm 


Figure  8.  Microfluidic  device  containing  BZ  droplets  connected  by  channels. 


We  also  succeeded  in  creating  so-called  double  emulsions,  which  offer  increased 
flexibility  in  linking  droplets  via  either  excitatory  or  inhibitory  coupling.  A  double  emulsion 
(W1/0/W2)  consists  of  an  outer  continuous  aqueous  phase  (wi),  a  shell  of  oil  (o),  and  an  inner 
aqueous  phase  (W2),  which  can  be  different  from  the  outer  phase.  In  our  case  the  two  aqueous 
phases  are  different  BZ  solutions  and  the  intermediate  oil  phase  is  a  fluorinated  oil.  The  process 
relies  on  proprietary  fluorinated  surfactants  synthesized  by  our  partner,  RDT.  Figure  9  shows 
double  emulsions  created  in  preliminary  experiments. 


Figure  9.  Double  emulsion 
(wi/o/w2).  BZ  solution  in  the 
oscillatory  state  is  in  the  core  of  the 
drop.  BZ  solution  in  the  excitable 
state  surrounds  the  drop.  Outer  shell 
of  drop  is  oil 
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3.  Theory 

The  systems  we  have  explored  basically  consist  of  arrays  of  coupled  oscillators,  where  each 
BZ  droplet  constitutes  a  chemical  oscillator,  and  the  coupling  arises  via  nonpolar  intermediates 
produced  in  the  course  of  the  reaction,  which  travel  through  the  oil  phase  between  adjacent 
droplets.  We  have  confirmed  by  addition  of  a  scavenging  molecule,  as  well  as  by  computer 
simulations,  that  the  species  that  acts  as  the  “messenger”  when  the  droplets  are  separated  by  oil 
is  elemental  bromine,  a  known  intermediate  of  the  BZ  reaction.  Coupled  oscillators  have  the 
potential  to  display  significantly  more  complex  behavior  than  individual  oscillators.  They  can 
generate  a  variety  of  phenomena,  the  simplest  of  which  involve  the  locking  of  phases  of  adjacent 
oscillators,  which  can  be  either  in-phase,  with  all  oscillators  acting  essentially  as  one,  or  out-of- 
phase,  typically,  but  not  necessarily,  by  180°,  so  that  adjacent  oscillators  are  alternately  “on”  or 
“off’.  If  one  can  understand  and  control  the  behavior  of  such  a  system  of  coupled  oscillators,  it 
should  be  possible  to  exploit  the  spatiotemporal  patterns  generated  not  only  for  chemical 
communication  but  potentially  for  chemical  computation  as  well. 

We  have  been  able  to  simulate  the  one-  and  two-dimensional  behavior  using  a  model  of  the 
BZ  reaction  that  takes  into  account  the  coupling  between  droplets  via  bromine.  The  simulations 
suggest  that  more  complex  patterns  can  arise.  Examples  of  ID  patterns  are  shown  in  Figure  10, 
and  Figure  1 1  shows  2D  patterns.  The  existence  of  multiple  attracting  states  suggests  that  it 
should  be  feasible  to  utilize  systems  of  this  type  for  complex  computational  as  well  as 
communications  applications. 
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Figure  10.  Exotic  regimes  found  in  simulations  of  6  coupled  oscillators  at  several  sets  of  model 
parameters.  For  all  simulations  initial  condition  is  ’’almost  in-phase,”  i.  e. ,  the  phase  of  one  of  the  six 
oscillators  is  slightly  shifted.  In  (a),  z\  and  z3  are  in-phase,  and  z4  and  z6  also  oscillate  in-phase.  In  (b),  z2 
and  z3  oscillate  in-phase,  and  z5  and  z6  are  also  in-phase,  while  zi  and  z4  are  almost  stationary  (small 
amplitude  oscillations).  In  (c)  z2  and  z5  oscillate  out-of-phase,  and  zi,  z3,  z4,  and  z6  are  almost  stationary. 
In  (d),  z2,  z3,  z5,  and  z6  behave  as  in  (b),  but  z\  and  z4  oscillate  in-phase  with  normal  amplitude  of 
oscillations  and  with  double  frequency.  zt  is  concentration  of  variable  corresponding  to  ferroin  in  the  i  ’th 
oscillator. 
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Figure  11.  Simulation  of  a  self-assembled  hexagonal  array  of  BZ-drops.  Left  panel  shows  a 
28-drop  array,  the  largest  we  have  simulated  to  date.  The  state  shown  is  the  “II-S”  state,  in 
which  one-third  of  the  drops,  shown  in  gray,  are  in  the  oxidized  steady  state,  one-third,  shown  in 
blue,  oscillate  with  nearly  the  same  phase,  and  the  other  one-third,  shown  in  red,  oscillated 
approximately  180°  out  of  phase  with  the  blue  drops.  Shades  of  blue  and  red  indicate  relative 
phases.  Differences  from  perfect  antiphase  behavior  result  from  differences  in  numbers  of 
nearest  neighbors  and  distance  from  the  boundary.  Right  panel  shows  time  series  of  three  blue 
drops,  1,  5  and  26,  illustrating  departure  from  perfect  synchrony. 

4.  Conclusion 

This  aspect  of  the  project  is  probably  the  most  promising  for  future  development.  The 
rich  behavior  of  arrays  of  coupled  oscillators  offers  a  wealth  of  possibilities,  both  scientific  and 
practical,  that  we  have  only  begun  to  explore.  While  implementation  in  communications  devices 
remains  a  challenge,  particularly  if  signals  must  be  transmitted  at  tens  to  hundreds  of  Hz, 
considerable  information  density  should  be  attainable.  More  appealing,  though,  are  applications 
to  “chemical  computing,”  where  it  should  be  possible  to  establish  the  “rules”  that  govern  how  the 
time  evolution  of  a  droplet  depends  on  the  states  of  its  neighbors,  i.e.,  one  could  construct  a 
“chemical  cellular  automaton.”  By  utilizing  the  photosensitivity  of  certain  versions  of  the  BZ 
system,  it  should  be  possible  to  build  and  control  an  analog  chemically-driven  computational 
device,  which  would  offer  advantages  for  certain  types  of  exponentially  difficult  (NP-complete) 
problems,  such  as  routing,  scheduling  and  data  retrieval.  We  have  developed  both  the 
experimental  and  the  theoretical  capabilities  to  understand  these  systems  in  considerable  detail. 

D.  Encoding 

We  explored  two  approaches  to  the  problem  of  encoding  the  signal  generated  by  the 
chemical  reactions.  We  report  our  results  below. 
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1.  Crack  and  Wrap 


The  scheme  we  denote  as  “crack  and  wrap”  involves  essentially  a  single  spatial  dimension. 
The  reactants  are  stored  in  separate  compartments,  and  the  reaction  is  initiated  by  breaking  a 
barrier  that  separates  them,  much  as  occurs  in  glowsticks  or  chemical  cold  packs.  The  device  is 
wrapped  in  a  field-printed  mask  that  modulates  the  visibility  of  the  periodic  chemistry  into  a 
coded  signal  from  which  the  message  can  be  recovered.  The  basic  scheme  is  shown  in  Figure 
12:  six  channels  are  used  to  encode  six  different  frequencies  in  order  to  give  the  six  bits  required 
to  encode  38  different  characters.  In  the  examples  shown  here,  the  chemiluminescent  solution  is 
a  commercial  glow  stick  and  the  mask  is  printed  on  a  transparent  mylar  sheet.  The  intensity  is 
measured  by  signal  averaging  and  Fourier  transformation  of  the  light  from  the  drops  propelled 
through  the  circuit.  In  Figure  13  the  intensity  of  light  from  one  drop  is  focused  on  the 
photomultiplier  detector  (PMT).  The  drops  generate  a  square  wave  signal  of  light  intensity  as  a 
function  of  time  as  shown  in  Fig.  13.  We  generated  signals  at  65  Hz  using  syringe  pumps  to 
propel  the  drops.  Figure  14  illustrates,  with  signals  from  two  channels,  our  ability  to  extract 
multiple  frequencies,  and  temporal  multiplexing  allows  us  to  send  sequences  of  characters.  The 
drops  can  be  transported  past  the  mask  either  by  an  external  power  source  or  by  chemical  waves 
(see  below). 


Figure  12.  (a)  Luminescent  drops  of  300  pm  diameter  in  cylindrical  tubing,  (b)  magnified  view  of  one 
bend,  (c)  a  single  drop  moving  under  a  mask  to  produce  a  periodic  modulation  of  the  light  intensity. 
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Figure  13.  (a)  Intensity  from  luminescent  drops  of  300  pm  diameter  moving  at  65  Hz  through  channels. 
The  detector  was  focused  on  one  drop,  (b)  Fourier  transform  of  light  intensity. 


Figure  14.  (a)  Different  masks  over 
each  of  the  six  lines  shown  in  Fig  4. 

(b)  Intensity  vs.  time  of  drops  moving 
under  one  mask  at  a  time  (green,  blue) 
and  under  two  masks  simultaneously 
(red),  (c)  FFT  of  red  trace  in  (b).  The 
FFT  shows  two  frequencies, 
corresponding  to  drops  simultaneously 
moving  under  two  masks. 
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2.  Spray  and  Display 

The  second  scheme,  “spray  and  display,”  utilizes  two  spatial  dimensions.  Here  the 
reactants  of  the  photosensitive  Ru(bpy)3-catalyzed  BZ  reaction  and  the  supporting  emulsion  are 
mixed  in  an  aerosol  can  containing  multiple  chambers  and  microfluidic  circuitry  and  propelled 
onto  a  prepared  surface,  where  they  self-assemble  into  a  regular  lattice,  e.g.,  hexagonal.  The 
technology  is  similar  to  that  needed  to  produce  the  foam  that  comes  out  of  a  shaving  cream  can. 
Alternatively,  the  target  surface  can  be  pre-pattemed  with  a  desired  array  of  drop  locations  and 
connections  as  in  Fig.  8.  With  the  photosensitive  reaction,  the  initial  condition  can  be  set  by 
projecting  the  desired  pattern  onto  the  lattice.  This  “chemical  cellular  automaton”  will  then 
evolve,  producing  a  long-lasting  repetitive  pattern  that  corresponds  to  the  message.  The  set  of 
initial  conditions  and  patterns  associated  with  the  array  of  possible  messages  can  be  developed 
through  extensive  simulations  and  experimentation  and  stored  in  a  codebook. 

3.  Conclusion 

The  “crack  and  wrap”  approach  appears  feasible  for  a  ID  system.  Considerable 
additional  work  needs  to  be  done  to  establish  the  feasibiity  of  the  “spray  and  display”  approach. 

E.  Power  Source 

In  order  for  approaches  involving  luminescent  materials  flowing  behind  masks  to  be 
successful,  it  is  necessary  to  have  a  means  of  creating  a  flow  of  fluid  through  the  system.  We 
have  explored  several  approaches  to  this  problem. 

1.  Gas  Pressure 

Our  initial  proposal  was  to  propel  material  through  the  microfluidic  system  by  generating 
a  gas  that  would  create  sufficient  pressure,  i.e.,  to  construct  a  device  analogous  either  to  an  air 
bag  or  a  pressure  cooker.  On  further  consideration,  we  concluded  that  such  approaches  are  likely 
to  be  subject  to  a  variety  of  technical  problems  involving  the  ability  to  control  the  gas  pressure. 

2.  Chemical  Waves 

A  much  simpler  approach  is  to  initiate  a  chemical  wave  that  would  carry  the  information, 
obviating  the  necessity  for  a  bulk  flow  of  material.  We  have  been  able  to  generate  waves  not 
only  in  continuous  media,  but  also,  as  illustrated  in  Fig.  6,  in  one-  and  two-dimensional  arrays  of 
droplets.  Whether  this  approach  is  compatible  with  the  desired  information  density  remains  to 
be  established. 

3.  Other  Approaches 

In  our  experiments,  we  used  syringe  pumps  to  drive  fluid  through  the  system.  This  approach 
is  clearly  not  adaptable  to  the  size  of  device  that  we  ultimately  sought  to  build.  Two  other 
approaches  do  seem  possible.  A  relatively  simple  solution  is  to  employ  printable  zinc-air 
microbatteries  (http://ccsl.mae.comell.edu/research/sff/index.htm  ),  which  can  provide  sufficient 
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power  at  very  little  cost  in  terms  of  added  mass  or  space  occupied.  It  may  also  be  possible  to 
incorporate  magnetic  nanoparticles  into  the  fluid  and  then  use  a  spinning  magnet  driven  by  a  tiny 
motor  to  generate  the  force  to  move  the  fluid. 

4.  Conclusion 

This  area  is  the  one  in  which  we  have  made  the  least  progress.  If  our  approaches  were  to 
be  carried  further,  considerably  more  work  would  need  to  be  devoted  to  power  sources.  The 
most  promising  approach  appears  to  be  the  use  of  small  batteries,  perhaps  of  the  zinc-air  type 
noted  above. 

F.  Signal  Detection 

The  SSI  and  AFRL  members  of  the  team  focused  their  efforts  on  the  issue  of  how  the 
signals  generated  by  the  device  can  be  detected  in  an  environment  that  contains  a  high  level 
of  temporally  varying  optical  noise.  To  deal  with  this  issue,  they  a)  developed  a  programmable 
simulator  to  mimic  the  signals  to  be  generated  by  the  chemical  communications  device;  and  b) 
carried  out  field  tests  to  assess  how  bright  a  modulated  source  might  have  to  be  to  enable  us  to 
distinguish  it  from  temporal  components  of  naturally  varying  backgrounds.  We  discuss  progress 
in  each  of  these  areas  below. 

1.  Programmable  simulator 

The  basic  design  of  the  simulator  is  shown  in  Figure  15.  It  consists  of  an  array  of  six 
infrared  LEDs.  Each  LED  is  individually  controllable  for  brightness  and  modulation,  and  there 
is  one  controller  channel  per  LED  and  one  analog  voltage  output  from  the  PC  for  each  controller 
channel.  Band  pass  filters  (BP)  are  used  to  select  wavelengths  not  at  the  LED  peak.  Neutral 
density  fdters  (ND)  provide  an  additional  level  of  power  control,  and  the  diffuser  (DF)  is 
available  to  modify  the  directional  LED  output. 
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2.  Field  Measurements 

On  June  10,  2008,  we  took  advantage  of  an  opportunity  to  utilize  AFRL  instrumentation  to 
monitor  a  modulated  light  source  in  a  natural  setting.  The  purpose  of  the  measurements  was  to 
develop  an  understanding  of  how  bright  a  modulated  source  might  have  to  be  to  be  able  to 
distinguish  it  from  temporal  components  of  naturally  varying  backgrounds. 

A  suite  of  instrumentation  was  located  in  the  Yankee  Building  atop  Mt.  Washington, 
including  radiometers  with  a  cutoff  filter  at  0.85,  and  band  pass  filters  at  1.25±0.075  and  1.6  ±0.1 
pm,  respectively.  In  addition,  an  infrared  imager  (Sensors  Unlimited,  Inc.  SU640)  fitted  with  a 
1.6  pm  filter  was  utilized.  The  radiometer  signals  were  sampled  at  20  kHz.  The  SU640  imager 
was  operated  at  a  120  HZ  frame  rate  in  256  x  256  pixel  mode.  The  radiometer  fields-of-view 
were  all  9.1°,  and  the  imager  FOV  was  49.1°.  Figure  16  shows  the  view  from  the  Yankee 
Building  atop  Mt.  Washington  looking  westward  toward  the  parking  lot  of  the  Cog  Railway 
Station,  where  the  light  source  was  operated.  The  range  from  the  instrument  suite  to  the  parking 
lot  was  4.16  km 
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Figure  16.  Clear  LOS  view  from  inside  the  Yankee  Building  atop  Mt.  Washington,  looking  westward 
toward  the  parking  lot  of  the  Cog  Railway  Station  (the  beige  pixels  at  the  center  of  the  red  circle). 

The  light  source  was  a  100  W  halogen  spot  light  with  a  9°  beam  width,  equipped  with  a 
Variac  voltage  control.  When  operated  at  full  power,  the  bulb  is  about  80%  efficient  with  a 
spectral  output  that  looks  like  a  3000  K  black  body.  Total  in-band,  in-beam  power  is  about  7-8 
W  in  each  of  the  1.2  and  1.6  pm  spectral  band  passes.  The  source  was  mounted  behind  a  20” 
five-blade  three-speed  box  fan  equipped  with  a  variable  voltage  control  to  adjust  the  speed.  The 
maximum  modulation  rate  for  this  setup  was  80  Hz. 

a.  Sample  results 

Figure  17  shows  the  raw  radiometer  signal  for  the  1.2  pm  band  pass.  For  this  particular 
experiment,  the  light  source  was  chopped  at  72  Hz  and  turned  on  full  power  at  t~10  sec,  then  the 
source  intensity  was  gradually  decreased  over  the  course  of  tens  of  seconds,  then  turned  on  back 
to  full  power  at  t~-T70  sec.  There  is  no  clear  indication  of  any  of  these  events  in  Figure  17.  Over 
the  course  of  this  run,  there  were  scattered  clouds  above  summit  altitude  passing  through  the 
scene.  All  of  the  temporal  variations  in  the  figure  are  due  to  time  varying  shadows  in  the 
radiometer  FOV.  Since  the  radiometer  FOV  is  quite  wide  (9.7°)  and  the  clouds  were  at  about 
summit  level,  cloud  movement  can  cause  gross  changes  in  the  total  signal  observed  in  the  space 
of  seconds.  Similar  behavior  was  observed  in  the  1.6  pm  channel  (Figure  18). 
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RAD  AC  06/10/2008  CH3  Signal 


Figure  17.  Raw  radiometer  signal  for  the  1 .2  pm  band  pass. 


RAD  AC  06/10/2008  CHI  Signal 


Figure  18.  Raw  radiometer  signal  for  the  1 .6  pm  band  pass. 
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In  order  to  pull  out  the  modulated  source  contribution  to  the  full  raw  radiometer  signal,  we 
typically  process  the  data  in  the  following  way.  First  the  data  is  de-trended  by  subtracting  off  the 
running  average  from  a  one  second  smoothing  window,  effectively  removing  the  DC  component 
as  well  as  frequencies  below  1  Hz  from  the  power  spectrum.  Then  a  Fourier  transform  is 
performed  on  the  data  in  1  sec  chunks.  The  result  is  a  power  spectrum  as  a  function  of  time  with 
low  frequencies  removed.  Figure  19  shows  the  result  using  the  radiometer  data  of  Figure  17.  In 
this  display,  fixed  frequency  events  show  up  as  vertical  lines.  In  Figure  19,  there  is  a  vertical  line 
at  72  Hz  that  shows  up  at  the  top  of  the  figure  and  decreases  in  intensity  with  time  until  it  is 
indistinguishable  from  the  background,  then  appears  again  at  170  s  when  the  source  is  returned 
to  full  power.  Figure  20  shows  similar  behavior  for  the  1.6  pm  channel.  The  difference  here  is 
that  there  appears  to  be  much  less  background  interference  at  lower  frequencies  at  1.6  pm  than 
for  1.2  pm,  with  the  result  that  the  source  at  72  Hz  is  easier  to  pick  out.  (The  prominent  feature  at 
60  Hz  is  due  to  pickup  of  line  voltage  at  some  point  in  the  instrument  detection  stream.) 


RAD  AC  06/10/2008  CHI  log  PSD 
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Figure  19.  Log  of  the  power  spectral  density  as  a  function  of  time  for  the  1 .2  pm  radiometer  channel. 
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Figure  20.  Log  of  the  power  spectral  density  as  a  function  of  time  for  the  1 .6  pm  radiometer  channel. 

Figure  21  shows  a  sample  frame  of  processed  data  from  the  SU640  imager.  Processing  of  this 
data  was  performed  as  for  the  radiometers,  except  that  it  was  performed  pixel  by  pixel  for  each 
pixel  in  the  display.  We  also  include  a  spatial  averaging  over  ~3  pixels  to  suppress  temporal 
changes  caused  by  pixel  to  pixel  movement  of  scene  features.  The  results  are  presented  as  a  quad 
with  the  display  at  the  upper  left  being  raw  imager  data,  and  the  display  at  the  upper  right  being 
an  RMS  image.  At  the  lower  left  is  pixel  power  summed  over  a  range  of  frequencies,  in  this  case, 
wide  enough  to  select  the  frequency  of  the  modulated  light  source.  Finally,  at  the  lower  right  is  a 
combination  of  vectors  resulting  from  principal  component  analysis  of  the  data.  Though  there  is 
no  indication  of  the  light  source  in  the  raw  image,  it  is  clearly  apparent  as  the  bright  yellow  spot 
just  over  the  rock  ledge  in  the  RMS  image.  The  source  is  also  visible  in  the  corresponding 
position  of  the  frequency  slice  at  the  lower  left,  though  it  is  not  so  obvious,  as  there  are  several 
interferences,  such  as  a  patch  of  sunlight  traveling  over  the  rock  ledge  in  the  near  field  at  the 
same  time.  The  PCA  image  at  the  lower  right  shows  little  or  no  evidence  of  the  source. 
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Imager  Data:  1213115538.9581  16:32:19 


RMS  Image:  1213115538.9581  16:32:19 
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Figure  21.  Quad  display  of  processed  output  from  the  SU640  infrared  imager  (see  text). 


3.  Conclusion 

The  simulator  provides  a  valuable  tool  for  testing  detection  capabilities,  and  the  ability  of  the 
SSI  team  to  carry  out  field  tests  and  analyze  signals  generated  should  prove  valuable  in  assessing 
the  capabilities  of  whatever  method  is  eventually  chosen  to  implement  the  desired  goals  for 
chemical  communication. 

III.  Publications 

This  work  has  resulted  in  three  publications: 

M.  Toiya,  V.K.  Vanag  and  I.R.  Epstein,  “Diffusively  Coupled  Chemical  Oscillators  in  a 
Microfluidic  Assembly,”  Angew.  Chem.  Int.  Ed.  47,  7753-7755  (2008) 

T.  Bansagi  Jr.,  M.  Leda,  M.  Toiya,  A.  M.  Zhabotinsky,  and  I.  R.  Epstein,  “High-frequency 
oscillations  in  the  Belousov-Zhabotinsky  reaction,”  J.  Phys.  Chem.  A  1 13,  5644-5648  (2009) 
M.  Toiya,  H.O.  Gonzalez-Ochoa,  V.  K.  Vanag,  S.  Fraden,  and  I.R.  Epstein,  “Synchronous 
Regimes  of  Oscillatory  Microdroplets,”  J.  Phys.  Chem.  Lett.  !,  1241-1246  (2010). 


Copies  of  these  works  appear  in  the  Appendix  below. 
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IV.  Appendix  -  Publications 
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Synchronized  Nano-Oscillators 


DOI :  1 0.1 002/anie. 20080233  9 


Diffusively  Coupled  Chemical  Oscillators  in  a  Microfluidic  Assembly** 

Mamhiro  Toiya ,  Vladimir  K.  Vanag ,  and  Irving  R.  Epstein * 


From  fireflies  that  synchronize  their  flashes  with  each  other11 1 
to  heart  muscles  contracting  and  relaxing  in  unison, p] 
synchronized  behavior  of  living  cells  or  organisms  is  ubiq¬ 
uitous  in  nature  J 31  Chemical  react  ton -diffusion  systems  can 
help  us  understand  the  mechanisms  that  underlie  such 
synchronization.  Coupled  chemical  oscillators  have  previ¬ 
ously  been  studied  in  the  laboratory  with  large  reactors 
connected  directly  by  small  channels  for  controlled  mass 
exchange  of  bulk  solution.*4-71  In  this  ease,  coupling  occurs  via 
all  species.  In  living  systems,  however,  coupling  often  occurs 
through  special  signaling  molecules,  as  in  synaptic  communi¬ 
cation  or  ehemotaxis.*81  Collections  of  neural  oscillators  can 
access  a  vast  repertoire  of  coordinated  behavior  by  utilizing  a 
variety  of  topologies  and  modes  of  coupling,  including  gap 
junctions  and  synaptic  links,  which  may  be  either  excitatory  or 
inhibitory,  depending  on  the  neuro  transmit  ter  involved. 

To  mimic  such  a  fine  level  of  communication  in  a  chemical 
system,  we  need  to  do  two  things:  a)  reduce  the  size  of  each 
oscillator  in  order  to  bring  the  characteristic  lime  of 
communication  between  diffusively  coupled  oscillators  to  or 
below  the  period  of  oscillation ;  and  b)  introduce  a  semi- 
permeahle  membrane  or  other  medium  between  tile  micro  - 
oscillators  to  permit  communication  only  via  selected  species. 
These  goals  can  be  achieved  with  the  use  of  microfluidic 
devices.  Our  experimental  system  (Figure  I  a)  is  a  linear  array 
of  tens  of  droplets  of  nanoliter  volume  containing  aqueous 
lerro  in-catalyzed  Belousov-Zhaboiinsky  ( BZ  solution  sep¬ 
arated  by  octane  drops  in  a  glass  capillary.  Ihe  BZ  reaction,  in 
which  the  oxidation  of  m  a  Ionic  acid  IMA)  by  hromale  is 
catalyzed  by  a  metal  complex  in  acidic  aqueous  solution,  is  a 
well  known  chemical  oscillator.  Owing  to  the  small  spatial 
extent  (/.  -  100-400  pm)  of  the  BZ  droplets,  the  characteristic 
lime  of  diffusive  mixing  within  a  single  droplet,  Iw2/D  (5- WO  s, 
D  —  diffusion  constant  of  aqueous  species),  is  smaller  than  the 
period  of  oscillation  ( 1 80-300  s) s  and  individual  BZ  droplets 
can  be  considered  homogeneous.  Bromine,  an  inhibitory 
intermediate  of  the  BZ  reaction,  is  quite  hydrophobic  and 
diffuses  readily  into  hydrocarbons  such  as  octane,  thus 
mediating  inhibitory  interdroplet  coupling.  We  have  shown 
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Figure  \  a)  Schematic  representation  of  the  microfluidic  device.  Red 
droplets  correspond  to  the  reduced  form  of  the  catalyst  {Ferroin),  blue 
droplets  to  the  oxidized  form  {ferriin).  A  new  method  for  fabricating 
such  junctions  is  outlined  in  the  Supporting  Information,  b)  Snapshot 
of  two  capillaries  with  droplets.  BZ  droplets  with  convex  surfaces  are 
dark  due  to  ferroin.  Horizontal  length  of  the  frame  and  inner  diameter 
{ID)  of  the  capillary  are  4.S  mm  and  150  pm.  respectively.  BZ  droplets 
were  recorded  by  a  CCD  camera  through  a  microscope  with  illumina¬ 
tion  by  light  passed  through  a  510  nm  interference  filter. 


theoretical  I  v1' 01  that  jn  such  heterogeneous  systems  patterns 
analogous  to  the  'luring  patternsJlM21  found  in  homogeneous 
systems  can  emerge. 

Without  comp  a  rime  ntaliza  lion,  the  homogeneous  BZ 
solution  in  a  similar  capillary  exhibits  trigger  waves  of 
excitation.  Partitioning  the  medium  into  droplets  dramati¬ 
cally  changes  this  behavior.  For  BZ  droplets  (Figure  1  b)  with 
/w>  400  pm  or  oil  droplets  with  length  /Q>400pm,  no 
discernible  coherent  patterns  are  seen.  However  if  ^*=100- 
400  pm  and  50-400  pm,  we  observe  stable  anti -phase 
oscillations  (Figure 2a)  at  larger  [MAf  (greater  than  100 mu) 
and  Turing  patterns  (Figure  2b)  at  smaller  [MA{  (less  than 
40  him).  At  higher  levels  of  [MA] ,  initially  in-phase  arrays  of 
droplets  evolve  to  an  anti  phase  configuration  within  a  few 
periods  of  oscillation  ( Movie  in  Supporting  Information).  For 
[M A]  — 40 him,  the  transition  to  the  Hiring  regime  goes 
through  intermediate  anti -phase  oscillations.  For  s  Eight  ly 
smaller  [MAJ  (35  him),  initially  in-phase  droplets  transform 
into  'lb ring  patterns  almost  immediately,  without  intermedi¬ 
ate  anti  phase  oscillations.  At  small  [MA[,  the  behavior  is 
rather  sensitive  to  the  size  of  droplets,  with  small  drops 
reaching  stationary  state  more  rapidly  than  larger  ones. 

To  establish  whether  bromine  is  responsible  for  commu¬ 
nication  between  the  BZ  droplets,  surfactant  Span  SO  (sorhi- 
tan  mo  no -o  I  eat  e)  at  concentrations  of  5  %  was  added  to  the 
octane.  In  separate  experiments,  it  was  found  that  Span 80, 
which  possesses  an  unsaturated  double  bond  in  its  hydro¬ 
carbon  tail  reacts  with  bromine  in  octane  in  less  than  1  s,  The 
water- insoluble  Span 80  thus  acts  as  a  trap  for  bromine, 
removing  il  from  the  octane.  When  Span 80  is  added  to  the 
droplet  system,  inhibiting  the  communication  between  BZ 
droplets,  individual  droplets  oscillate  independently,  if  we 
initiate  the  system  (see  Experimental  Section)  with  all 
droplets  in  the  same  phase,  in-phase  oscillations  persist. 
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Figure  2.  Space-time  plots  showing  a)  anti-phase  oscillations  with 
spikes  of  oxidation  offerroin  seen  as  light  horizontal  lines  across  BZ- 
droplets  and  b)  stationary  Turing  structures  with  alternating  oxidized 
and  reduced  states  evolving  from  an  initial  oscillatory  state.  Horizontal 
lengths  of  the  frame  and  the  capillary  ID  are  4.8  mm  and  150  Jim, 
respectively;  the  total  times  for  (a)  and  (b)  are  5200  s  and  10800  s, 
respectively.  Patterns  extend  to  the  left  and  right  of  the  segments 
shown. 


Seeking  to  understand  our  experimental  results  in  more 
depth,  we  performed  a  series  of  computer  simulations  to 
ascertain  whether  bromine  is  the  single,  or  at  least  the  major, 
signaling  molecule,  and  whether  other  patterns  may  be  found 
at  smaller  droplet  sizes  inaccessible  in  our  experiments.  The 
BZ  reaction  generates  a  second,  excitatory  intermediate, 
Br02\  which  is  also  capable  of  diffusing  in  the  oil  phase.  The 
model,  which  is  described  in  detail  in  the  Supporting 
Information,  is  based  on  the  Field- Koros-Noyes  (FKN) 
mechanism^13!  of  the  BZ  reaction  and  employs  seven  concen¬ 
tration  variables  to  describe  the  aqueous  phase  and  two  more, 
corresponding  to  [Br2]  and  [BrOy],  for  the  oil  phase.  The 
concentrations  of  the  major  reactants,  H+,  MA,  BrMA,  and 
Br03"  which  are  significantly  larger  than  those  of  the  variable 
species,  are  taken  to  be  fixed  in  a  given  experiment.  The 
model  contains  9n  variables  for  n  coupled  oscillators.  In 
addition  to  the  initial  concentrations,  key  parameters  in  the 
model  include  the  coupling  constants,  k(  and  kb ,  which 
characterize  the  strength  of  coupling  mediated  by  Br2  and 
BrOy,  respectively.  We  simulated  arrays  with  two  (Figure  3a), 
four,  and  six  coupled  oscillators  to  investigate  how  the 
behavior  of  the  system  depends  on  the  number  of  oscillators. 

The  two  stable  modes  found  in  the  experiments,  anti¬ 
phase  oscillations  and  Turing  patterns,  are  seen  in  the  model 
at  large  and  small  [MA],  respectively,  and  are  shown  in 
Figure  3b  and  c.  Note  that  if  kf  =  0,  that  is,  coupling  via  Br2  is 
absent,  neither  the  anti-phase  oscillations  nor  the  Turing 
mode  occurs,  so  Br2  is  an  essential  “messenger"  for  these  two 
regimes.  At  higher  [MA]  (greater  than  0.2  m),  where  the  anti¬ 
phase  mode  dominates,  the  results  of  the  simulation  are 


o 


© 


[MA]/m 


Figure  3.  a)  Configuration  of  two  coupled  BZ  oscillators  used  in 
simulations,  b)  Typical  anti-phase  oscillations,  c)  Typical  Turing  mode; 
bold  and  dashed  lines  in  (b)  and  (c)  represent  z  in  two  neighboring 
BZ  droplets,  d)  Parametric  diagram  in  [MA]-kf  plane  for  two  oscillators 
whose  initial  phases  are  slightly  shifted  (kb=kf/Pu,  kbt=kb/Ps;  partition 
coefficients  Pu  =  20,  Ps  — 1,  [H *]  =  /?  =  0.2  m,  [Br03_]  =  A  =  0.3  m, 
z  +  c  =  3  mM) ;  D=2xKTs  cm2s-1  is  used  for  conversion  ofkf  into  a 
characteristic  length  L=  (2D/kf)1/2.  Symbols:  o  and  A:  Turing  mode  for 
kt  =  0  and  for  both  kfr  =  0  and  kt  =  kf,  respectively;  ■;  anti-phase 
oscillations;  unstable  in-phase  oscillations  transforming  into  anti¬ 
phase  oscillations  for  larger  number  (four  or  six)  of  coupled  oscilla¬ 
tors;  x:  stable  in-phase  oscillations  (for  kfr=0);  Turing  mode  marked 
by  o  is  replaced  by  in-phase  oscillations  (x)  for  kfl=k(;  +:  weak 
communication  (initial  phases  of  oscillators  change  very  slightly  after 
ten  periods  of  oscillations).  Subscripts  “12”  and  “21”  refer  to  U  ([Br2]) 
and  S  ([BrO/])  in  octane  droplets. 


nearly  independent  of  the  presence  of  Br02\  while  at  lower 
[M  A]  (less  than  0.1m)  and  at  large  kf  (greater  than  0.5  s-1),  the 
Turing  mode  dominates  at  kh  =  0  and  the  in-phase  mode 
dominates  at  kh  =  kf  (Figure  3d).  Since  at  kh  =  kf  the  Turing 
mode  is  found  only  at  [MA]<20mM,  while  in  our  experi¬ 
ments  we  found  the  Turing  mode  at  larger  [MA]  (40  m.vi),  we 
infer  that  Br02‘  plays  a  minor,  if  any,  role  in  communication 
between  the  BZ  droplets. 

For  many  sets  of  parameter  values,  two  or  more  modes  are 
simultaneously  stable,  and  the  mode  obtained  depends  upon 
the  initial  conditions.  Simulations  with  four  and  six  oscillators, 
however,  reveal  that,  when  they  coexist,  the  in-phase  mode  is 
always  less  stable  than  the  Turing  or  anti-phase  modes.  Stable 
in-phase  behavior  is  found  only  at  large  kh  —  k{,  corresponding 
to  very  small  droplet  lengths  (less  than  100  pm).  With  such 
small  droplets,  we  also  find  several  more  exotic  regimes,  some 
of  which  are  illustrated  in  the  Supporting  Information. 

Chemical  nano -oscillators  diffusively  coupled  by  known 
signaling  species  may  provide  useful  analogs  for  biological 
processes.  The  microfluidic  BZ-octane  system  employed  here 
is  convenient  in  that  we  are  able  to  identify  the  inhibitor 
bromine  as  the  main  messenger  species,  and  the  production 
and  function  of  bromine  in  the  overall  BZ  process  are  well 
characterized.  By  choosing  the  fundamental  oscillator  and  the 
scavenger  species  added  to  the  connecting  medium,  it  should 
be  possible  to  build  systems  with  controllable  degrees  of 
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inhibitory  or  excitatory  coupling.  Microfluidic  technology 
makes  it  possible  to  construct  two-  as  well  as  three-dimen¬ 
sional  arrays  of  coupled  oscillators.  The  possibility  of  devel¬ 
oping  computational  devices  by  combining  oscillatory  chem¬ 
ical  reactions  with  droplet-based  microfluidic  techniques  has 
recently  been  suggested.1141 

Experimental  Section 

BZ  mixture:  The  aqueous  reaction  mixture  contains  H2S04  (80  m.vi), 
NaBr03  (0.288vi),  and  ferroin  (3  mvi).  In  addition,  for  experiments  on 
the  oscillatory  mode  we  add  NaBr  ( 10  mvi)  and  MA  (0.64m),  while  for 
experiments  on  Turing  patterns  [NaBr]=0,  [MA]  =  30-50  mvi.  To 
make  the  BZ  reaction  photosensitive,  we  add  a  small  amount 
(0.4  mw)  of  [Ru(bpy)3]  (bpy  =  bipyridinc).(151 

Microfluidics:  The  BZ  solution  and  octane  are  driven  by  syringe 
pumps  into  a  microfluidic  junction  at  the  entrance  to  the  capillary  as 
shown  in  Figure  1  a.  The  sizes  of  and  separation  between  the  BZ 
droplets  depend  upon  the  junction  size  and  the  flow  rates  with  which 
the  two  components  arc  injected  into  the  systcm.fl61  Before  the  start  of 
each  experiment,  we  put  the  BZ  micro-oscillators  into  the  in-phase 
mode  by  illuminating  the  capillary  with  a  strong  450  nm  light.  Further 
technical  details  arc  given  in  the  Supporting  Information. 
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Chemical  oscillations  in  the  classic  Belousov— Zhabotinsky  (BZ)  system  typically  have  a  period  of  a  few 
minutes,  which  can  be  increased  significantly  by  changing  the  organic  substrate.  Here  we  show  that  by  changing 
the  temperature  and  concentrations,  an  increase  of  3—4  orders  of  magnitude  in  the  frequency  of  BZ  oscillations 
can  be  obtained.  At  elevated  temperatures,  in  high  concentration  mixtures,  the  cerium-catalyzed  reaction 
exhibits  sinusoidal  oscillations  with  frequencies  of  10  Hz  or  greater.  We  report  the  effect  of  temperature  on 
the  frequency  and  shape  of  oscillations  in  experiments  under  batch  conditions  and  in  a  four-variable  model. 
We  show  that  our  simple  model  accurately  captures  the  complex  temporal  behavior  of  the  system  and  suggests 
paths  toward  even  higher  frequencies. 


Introduction 

The  Belousov— Zhabotinsky  (BZ)  reaction1*2  continues  to 
draw  considerable  scientific  interest  even  five  decades  after  its 
discovery.  Over  the  years,  the  name  has  come  to  refer  to  a  family 
of  metal-ion-catalyzed,  bromate-driven  reactions  in  which 
organic  substrates  in  acidic  media  are  oxidized  to  give  oscil¬ 
lations  in  the  concentrations  of  various  intermediates.  The 
detailed  mechanism  of  the  classical,  cerium-ion-catalyzed  reac¬ 
tion  system  was  first  described  by  Field.  Koros,  and  Noyes3 
(FKN).  By  reducing  its  complexity  to  a  few  coupled  elementary 
pseudoreactions.  Field  and  Noyes  were  able  to  construct  a 
simplified  mathematical  model,4  now  known  as  the  Oregonator, 
which  exhibits  oscillatory  limit  cycles.  The  model  also  suc¬ 
cessfully  describes  a  wealth  of  complex  phenomena  seen  in  the 
BZ  reaction,  both  temporal,  such  as  excitability,5*6  bistability,7 
stirring  effects,8*9  quasiperiodicity,  and  chaos,10*11  and  spatial, 
including  traveling  waves,12  target  and  spiral  patterns,13*14  scroll 
rings,15  and  their  orientation  dynamics  in  advective  fields.16 

The  Oregonator  has  been  also  used  to  model  the  effect  of 
temperature  on  oscillations,  prompted  by  Koros’  experimental 
study17  of  the  BZ  reaction  between  15  and  35  °C  with  Ce3+, 
Mn2+,  and  Ru(bipy)32+  as  catalysts.  He  pointed  out  that  the 
complex  BZ  system  behaves,  in  a  sense,  as  if  there  is  a  single 
rate-determining  step,  since  it  obeys  the  Arrhenius  equation, 
yielding  practically  the  same  activation  energy  for  all  three 
catalysts.  In  addition  to  experimentally  confirming  these  striking 
findings,  Blandamer  et  al.,18*19  by  systematically  varying  the  rate 
constants  in  the  model,  identified  the  model  step  that  most 
influences  the  frequency  as  the  one  representing  the  reaction 
between  bromide  and  bromate  ions  to  yield  bromous  acid. 
Misra20  studied  the  BZ  reaction  both  experimentally  and 
numerically  over  a  significantly  wider  temperature  regime,  from 
30  to  90  °C,  in  a  continuous-flow  stirred  tank  reactor  (CSTR  ). 
By  using  an  extended  Oregonator  model  that  incorporates  the 
bromination  of  malonic  acid  by  hypobromous  acid  and  bromine 
and  the  oxidation  of  malonic  acid  and  its  brominated  derivatives 
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by  Ce4+  and  other  species  to  form  bromide  and  free  radicals, 
he  demonstrated  that  the  frequency  cannot  continue  to  grow 
exponentially  without  limit.  Instead,  with  increasing  temperature, 
oscillations  disappear  at  a  critical  temperature  via  a  Hopf 
bifurcation.  Over  this  same  range  of  temperatures,  the  amplitude 
of  oscillations,  however,  remains  essentially  unchanged  before 
quickly  dropping  near  the  bifurcation  point.  Similar  conclusions 
were  drawn  by  Strizhak  and  Didenko21  for  temperatures  below 
33  °C,  employing  the  kinetic  scheme  of  Gyorgyi  and  Field.22 

None  of  the  studies  above,  regardless  of  the  model  or  scheme 
employed,  examined  the  effects  of  varying  the  acidity  or  howr 
the  shape  of  the  oscillations  varies  with  temperature.  Concentra¬ 
tions  of  reactants,  except  for  sulfuric  acid,  were  typically 
maintained  in  the  range  explored  in  the  prior  literature  for  the 
BZ  reaction.  Although  increases  in  temperature  did,  as  expected, 
shorten  the  period  from  typical  values  of  ~5  min  to  as  low  as 
2  s,20  one  may  wonder  what  range  of  frequencies  is  accessible 
to  the  BZ  reaction.  Use  of  alkylated  species,23*24  such  as 
methylmalonic  acid,  or  aromatic  species,25  such  as  pyrocatechol. 
as  the  substrate  can  lengthen  the  period  by  a  factor  of  10  or 
more.  Might  it  be  possible  to  decrease  the  period  by  several 
orders  of  magnitude,  and  could  one  model  such  high-frequency 
oscillations  within  the  Oregonator  framework?  Here,  we  seek 
to  answer  these  questions  by  exploring  the  uncharted  regime 
of  higher  concentrations  and  temperatures,  aiming  to  find  high- 
frequency  oscillations  and  to  numerically  reproduce,  with  a  high 
degree  of  accuracy,  the  temperature  dependence  of  the  frequency 
and  amplitude  as  well  as  the  shape  of  oscillations  with  a  four- 
variable  Oregonator  model.  We  also  use  the  model  to  roam  the 
reaction’s  phase  space  in  search  of  even  higher  frequency 
oscillations  in  this  well-studied  chemical  system. 

Experimental  Section 

Oscillatory  profiles  are  obtained  by  electrochemical  measure¬ 
ments  within  the  temperature  range  of  40—80  °C.  Reaction 
mixtures  are  prepared  from  aqueous  stock  solutions  of  analyti¬ 
cal-grade  sodium  bromate  (NaBr03,  Fisher),  malonic  acid 
(CH2(COOH)2,  Sigma-Aldrich),  sulfuric  acid  (H2S04,  Fisher, 
5  M),  and  ceric  ammonium  nitrate  ((NH4)2Ce(N03)6,  Fisher). 
The  initial  concentrations  are  [NaBr03]  =  0.5  M,  [CH2- 
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Figure  1.  Oscillatory  traces  at  three  temperatures  in  the  BZ  reaction: 
40  (dotted),  60  (dashed),  and  80  °C  (solid).  Initial  conditions  are  given 
in  the  Experimental  Section.  Profiles  are  shifted  along  the  abscissa  for 
clarity. 


Figure  2.  Initial  oscillatory  frequencies  (squares)  and  amplitudes 
(triangles)  in  the  temperature  range  40—80  °C.  Dashed  line  represents 
exponential  curve  fitted  to  experimental  frequency  data.  Dotted  line 
corresponds  to  the  cubic  curve  drawn  to  aid  the  eye.  Reactant 
concentrations  are  listed  in  the  Experimental  Section. 


(COOH)2]  =  2.0  M,  [H2S04]  =  1.5  M,  [(N^hCefNCbk]  = 
1 .0  mM,  except  where  otherwise  specified. 

As  the  working  electrode  we  use  a  glass-coated  platinum  wire 
with  one  end  open  to  the  electrolyte,  immersed  in  0.25  mL  of 
BZ  mixture  held  in  a  thin-walled  5  mm  diameter  glass  tube. 
The  working  electrode  potential  is  measured  against  a  saturated 
calomel  electrode  (SCE,  Fisher  Scientific  Accumet,  gel  filled) 
placed  in  a  1 .5  cm  diameter  glass  tube  containing  0.25  M  sodium 
sulfate  (Na2SC>4,  Fisher)  solution.  Both  half  cells  are  submerged 
in  a  thermostat  (Techne  equipped  with  Techne  TE-10D  digital 
immersion  circulator)  and  connected  through  a  salt  bridge,  a  1 
mm  diameter  glass  U-tube  filled  with  0.25  M  sodium  sulfate  in 
an  agar— gel  matrix  (melting  point  %  85  °C).  Reaction  mixtures 
are  stirred  with  a  0.8  mm  diameter  glass  rod  attached  to  a  high- 
rpm  drill  (Dremel  300).  The  potential  difference  is  recorded 
with  a  high-performance  USB  acquisition  module  (Data  Trans¬ 
lation  9812)  at  sampling  frequencies  ranging  from  100  to  1000 
Hz.  Despite  the  significant  heat  production,  working  with  such 
small  volumes  of  reaction  mixture  in  the  arrangement  described 
above  enables  us  to  reliably  keep  the  temperature  within  a  ±1 
°C  range  around  its  target  value  in  the  BZ  half  cell.  Evaporation 
is  negligible,  as  there  is  no  discernible  change  in  the  height  of 
the  menisci,  even  at  80  °C.  We  are  able  to  resolve  oscillations 
with  amplitudes  of  1 0  mV  and  higher. 

After  lowering  the  fully  assembled  electrochemical  cell  into 
the  water  bath  of  the  thermostat,  we  inject  the  requisite  amounts 
of  all  preconditioned  reagents,  except  the  catalyst,  into  the  BZ- 
batch  half  cell.  Subsequently,  under  vigorous  stirring,  we  add 
the  catalyst  to  start  the  reaction. 

Experimental  Results 

The  preoscillatory  induction  period  decreases  with  temper¬ 
ature,  from  about  300  s  at  40  °C  to  about  30  s  at  80  °C.  During 
this  portion  of  the  reaction,  the  potential  gradually  decreases 
from  its  initial  maximum  value.  This  downward  drift  of  the 
average  potential  continues  throughout  the  oscillations  and 
beyond,  until  the  half-cell  potential  difference  reaches  zero.  This 
effect  can,  however,  be  neglected  within  sufficiently  short 
intervals. 

Typical  oscillatory  profiles  at  several  temperatures  are 
shown  in  Figure  1.  For  better  comparison,  each  is  displaced 
along  the  time  axis  to  fit  within  a  1  s  interval.  The  profiles 
reveal  the  strong  temperature  dependence  of  the  oscillatory 
frequency.  The  waveform  is  nearly  sinusoidal  between  40 
and  80  °C,  in  contrast  to  the  relaxation  oscillations  observed 
closer  to  room  temperature.  The  data  also  show  a  small 
decrease  in  amplitude  over  the  given  temperature  range. 


In  all  cases,  as  the  reaction  proceeds,  the  frequency  gradually 
decreases  until  oscillations  finally  cease  with  a  period  20—30% 
longer  than  at  the  start.  The  number  of  cycles  lessens  with 
increasing  temperature  and  drops  to  about  15  at  80  °C.  The 
relative  position  of  the  curves  along  the  ordinate  suggests  that 
the  average  potential  rises  with  increasing  temperature,  though 
this  conclusion  is  somewhat  uncertain  owing  to  small  variations 
( <50  mV)  from  experiment  to  experiment  in  the  overall 
potential  differences  between  the  BZ  and  SCE  half  cells. 

Initial  oscillatory  frequencies  and  amplitudes  for  temperatures 
between  40  and  80  °C  are  presented  in  Figure  2.  The  frequency 
grows  exponentially  with  temperature  over  the  studied  range, 
while  the  amplitude  peaks  at  about  65  °C.  The  calculated 
activation  energy  is  58.0  ±  0.7  kJ/mol,  in  good  agreement  with 
previous  studies. 17-19-26 

We  also  carried  out  experiments  w  ith  1,4-cyclohexanedione 
(CHD)  as  the  organic  substrate.  In  this  set  of  experiments, 
ferroin  is  used  as  the  catalyst  and  the  traces  are  acquired 
spectrophotometrically.  The  maximum  frequency  of  sinusoidal 
oscillations  is  measured  to  be  12  Hz  with  [CHD]0  =  0.6  M. 
[NaBrO3]0  =  0.4  M.  [H2SO4]0  =  2.0  M.  and  [Ferroin]  =  0.5 
mM  at  85  °C.  We  note  that  under  such  conditions  we  see 
surprisingly  vigorous  gas  production  in  this  reaction,27  which 
is  often  regarded  as  a  bubble-free  BZ  variant. 

Model 

The  numerical  analysis  is  based  on  a  model  incorporating 
all  reactant  species  proposed  by  Zhabotinsky  et  al.28  to  describe 
the  shape  of  oscillations  and  wave  properties  in  the  BZ  reaction. 
The  simplified  set  of  equations  used  to  describe  our  experimental 
findings  is 

dX/dt  =  —k1h0XY  +  k^Y- 2k/X2  -  k^fiX  +  k_5U 2  + 

k6u(c-Z)-k_ja  (i) 

d  Y/dt  =  —k2hgXY  -  kjtfi  Y  +  qk7BZ  +  k9B  (2) 

dZ/dt  —  k6U(C  —  Z)  —  k_6XZ  —  k7BZ  (3) 

dU/dl  =  IkfhgAX  -  2k_5U2  -k6U(C-Z)  +  k_JCZ  (4) 

where  X  =  |HBr02],  Y=  [Br "],Z  =  [Ce4+],  U  =  |HBr02]. //<, 
is  Hanunet’s  acidity  function,29  A  =  [NaBrC)3]  =  ho/(0.2  +  h0 ) 
[NaBrO3]0,  B  =  [MA],  C  =  [Ce3+]  +  [Ce4+],  and  k4'  =  k4(  1  + 
0.87/fo).  For  numerical  calculations,  based  on  the  initial  reactant 
concentrations,  we  assign  Ao  =  0.5  M.  Bo  =  2  M,  C  =  0.001 
M.  and  ho  =  3,  which  corresponds  to  [H2S04]  =  1.5  M.  The 
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TABLE  1:  Rate  Parameters® 


Bansagi  et  al. 


rate  constant 

20  °C 

40  °C 

AC  [kJ/mol] 

A* 

80  °C 

80  °C  (this  paper) 

At2  [M-2  s_l] 

7.57  x  106 

1.2  x  107 

15.0 

1.24  x  107 

2.60  x  107 

4.95  x  107 

k3  [M“l  s_I] 

2.0 

15 

74.3 

1.20  x  10* 1 

430 

81.7 

A'4  [M~2s_1] 

8.6  x  103 

1  x  104 

3.23 

111 

1.3  x  104 

2.5  x  104 

*5[i srV] 

10 

60 

65.8 

1.83  x  10‘° 

1187 

3158 

k-s  [M_l  s-1] 

4.2  x  106 

6  x  106 

11.1 

1.35  x  106 

1.1  x  107 

2.1  x  107 

MM-'s”1] 

5  x  104 

9  x  10* 

19.9 

6.00  x  105 

2.4  x  10s 

4.6  x  10 5 

A'-6  [M-1  s"1] 

5  x  104 

2  x  104 

5.99 

798 

3.7  x  104 

6.95  x  104 

k7  [M-1  s_l] 

0.4 

0.8 

23.9 

25.0 

2.56 

30.5 

Ms"1] 

3.3  x  10“6 

5  x  10"6 

13.3 

2.67  x  10"6 

1.0  x  10"5 

2.4  x  10“5 

a  Values  in  columns  2  and  3  are  taken  from  ref  28.  Values  of  the  energy  of  activation  AG  in  column  4  and  the  pre-exponential  factor  A  in 
column  5  are  computed  from  values  of  the  rate  constants  in  columns  2  and  3.  Values  in  column  6  are  obtained  by  extrapolation  using  values 
from  columns  2  and  3.  In  column  7,  values  are  adjusted  to  experiments  at  80  °C.  b  Units  are  same  as  those  for  kt. 


Figure  3.  (a)  Oscillatory  solutions  to  systems  1—4  at  T  —  80  (black). 
60  (gray),  and  40  °C  (dotted).  A0  =  0.5  M,  B  =  2  M,  /i0  =  3  M,  and 
C  =  0.001  M.  (b)  Projections  of  limit  cycles  (composed  of  circles)  on 
phase  plane  log  X  vs  ZJC  for  temperatures  T  =  80  (blue),  60  (green), 
and  40  °C  (red).  Time  elapsed  between  consecutive  symbols:  5  ms. 
Solid  and  dashed  lines  are  zero-growth  isoclines  (nullclines)  for  X  and 
Z.  respectively,  calculated  for  the  three  temperatures.  Note  that  the 
abscissa  is  logarithmic. 

stoichiometric  parameter  q  is  set  to  0.6  in  accordance  with  ref 
28. 

Numerical  Results 

Modeling  our  experimental  findings  requires  rate  constants 
applicable  al  high  temperatures.  For  the  reactions  that  lead  to 
eqs  1—4.  however,  these  quantities  are  available28  only  for  20 
and  40  °C  (Table  1,  columns  2  mid  3).  To  estimate  their  values 
at  higher  temperatures,  we  can  use  the  absolute  rate  theory, 
according  to  the  which  a  rate  constant  can  be  calculated  as  lc(T) 
=  AT exp(— EJRT).  where  k.  A.  T.  and  £a  are  the  rate  constant, 
pre-exponential  factor,  absolute  temperature,  and  activation 
energy,  respectively.  For  the  rate  constant  values  given  by  this 
formula  (listed  in  column  6  of  Table  1),  oscillations  vanish  at 
temperatures  as  low  as  55  °C  with  frequencies  below  1  Hz  at 
the  Hopf  bifurcation.  In  our  experiments,  however,  we  observe 
oscillations  at  60  °C  (~3  Hz)  and  higher  temperatures. 

To  address  this  discrepancy  and  to  quantitatively  capture  the 
experimentally  observed  oscillations,  we  adjusted  all  rate 
constants  to  obtain  the  best  fit  with  the  experimental  data, 
specifically  the  frequency  and  amplitude  of  oscillations  at  80 
°C  for  A0  =  0.5  MJ0  =  2M.C=  0.001  M,  and  hG  =  3  and 
the  concentration  of  bromate  A0  at  which  the  Hopf  bifurcations 
at  40  and  60  °C  occur  for  Bo  =  2  M.  C  =  0.001  M.  and  ho  = 
3  (Figure  5a).  We  chose  the  bromate  vs  temperature  phase 
diagram  (Figure  5a)  because  the  system  is  quite  sensitive  to 


40  50  60  70  80  90 

temperature  (°C) 


Figure  4.  Simulated  dependences  of  frequency  (squares)  and  amplitude 
(triangles)  of  potential  oscillations  on  temperature.  Amplitude  is 
calculated  as  /?77F[ln(Zmax/(C  -  Z^,))  -  ln( Z^/tC  -  Zaun))],  where 
F  is  Faraday’s  constant.  Parameter  values  are  the  same  as  in  Figure  3. 


changes  in  the  bromate  concentration.  We  believe  that  this 
approach  is  more  reliable  than  simply  lilting  potential  time  series 
because  the  fast  oscillations  last  only  a  few  seconds. 

Comparing  the  resulting  constants  (column  7)  to  their 
counterparts  given  by  the  absolute  rate  theory  (column  6),  we 
find  that  only  the  differences  in  the  k?,  and  k-j  values  (rows  3 
and  9)  are  significant.  These  steps  correspond,  respectively,  to 
the  following  reactions 

H+  4-  Br~  +  HBrO,-*-’-  HBi02  +  HOBr  (5) 
Ce4+  +  CHBr(COOH),— Ce3+  +  CBitCOOH)j  +  H+ 

(6) 

Reaction  5  was  identified  by  Blandamer  at  al.19  as  the  major 
determinant  of  the  oscillatory  frequency,  and  its  activation 
energy  was  measured  to  be  54  kJ/mol.30  Substituting  this  value 
and  the  corresponding  rate  constant  at  20  °C  (Column  2)  into 
the  absolute  rate  theory  equation  yields  a  rate  constant  of  86.6 
M-1  s"1  at  80  °C.  which  is  in  excellent  agreement  with  that 
calculated  in  this  study.  The  activation  energy  of  reaction  6  is 
reported  to  be  67  kJ/mol21  and  is  assumed  to  be  close  to  the 
total  energy  of  activation,31  which  is  58  kJ/mol.  With  the 
extrapolation  method,  these  two  values  of  E.d  both  yield  fairly 
good  estimates  for  kj .  42.9  and  22.9  M-1  s-1,  respectively. 
Hence,  we  conclude  that  our  corrections  for  constants  k^  and  A'7 
are  reasonable  and  consistent  with  literature  data. 

With  these  adjustments  to  the  rate  constants,  the  calculated 
frequencies,  amplitudes,  and  shapes  of  the  oscillatory  profiles 
faithfully  reflect  our  experimental  findings  (see  Figure  3a). 
Analyzing  the  dynamics  of  the  four- variable  model  on  the  phase 
plane  provides  detailed  insights  into  the  behavior  of  the  BZ 
reaction  at  elevated  temperatures.  According  to  Figure  3b,  the 
limit  cycle  is  fairly  symmetric  for  all  three  temperatures,  which 
is  consistent  with  the  sinusoidal  character  of  the  oscillatory 
traces.  As  the  temperature  is  increased,  the  limit  cycle  shrinks 
and  the  oscillations  become  substantially  faster.  The  distance 
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Figure  5.  Calculated  parametric  diagrams  T—Ao  (a),  T—B  (b),  and 
T—h0  (c).  Remaining  parameters  are  the  same  as  in  Figure  1.  Black 
lines  indicate  boundaries  between  dynamic  states.  Symbols  OSC  and 
SS  denote  areas  corresponding  to  oscillation  and  stationary  states, 
respectively.  The  secondary  ordinate  axis  shows  the  oscillation 
frequency  (gray  lines)  at  Hopf  bifurcation  points.  Circles  correspond 
to  experimental  results:  no  oscillations  (empty  circles),  oscillations  (solid 
back  circles),  with  frequencies  shown  by  solid  gray  circles. 

between  the  extrema  of  Z,  however,  changes  much  less  rapidly, 
which  translates  to  oscillations  with  nearly  constant  amplitude 
throughout  the  studied  temperature  regime,  as  observed  in 
experiment  and  previously  reported  for  dilute  BZ  mixtures.20 

The  oscillatory  frequency  and  amplitude  dependence  on 
temperature  obtained  from  the  model  are  presented  in  Figure 
4.  As  in  the  experiments,  there  is  an  exponential  relationship 
between  frequency  and  temperature.  The  calculated  amplitude 
exhibits  the  same  tendencies  seen  in  the  experiments  (Figure 
2).  Compared  to  the  frequency,  the  amplitude  remains  within  a 
narrow  range,  lirst  grow  ing,  then  passing  through  a  maximum, 
and  finally  decreasing  as  the  temperature  approaches  80  °C. 
Simulations  suggest  that  the  sharp  fall  in  amplitude  continues 
while  the  frequency  rapidly  increases  until  oscillations  finally 
disappear  around  90  °C  with  a  frequency  of  about  20  Hz. 
Spectrophotometric  measurements  at  higher  temperatures  with 
ferroin  as  catalyst  using  slightly  different  reactant  concentrations 
support  this  prediction.  The  oscillator)'  regime  in  these  ferroin- 
catalyzed  experiments  extends  above  80  °C,  and  the  frequency 
can  reach  as  high  as  18  Hz  for  temperatures  near  100  °C. 

The  fitted  total  activation  energy  is  68.2  ±  0. 1  kJ/mol,  which 
exceeds  that  found  in  experiment.  This  difference  is  a  conse¬ 
quence  of  the  tendency  of  the  model  to  underestimate  frequen¬ 
cies  for  lower  temperatures,  as  it  is  primarily  designed  to 
accurately  describe  the  high-frequency  oscillations. 

Our  model,  despite  the  aforementioned  weakness,  enables  us 
to  numerically  explore  the  high-frequency  oscillatory  regime 
of  the  BZ  system.  When  mapping  this  domain  of  the  reaction’s 
phase  space,  we  keep  all  but  one  parameter  as  well  as  the  total 
catalyst  concentration  constant  while  systematically  changing 
the  remaining  parameter  and  the  temperature.  Guided  by  the 


numerical  results,  we  carry  out  experiments  along  sections  of 
the  calculated  boundaries  separating  oscillatory  and  stationary 
regimes  to  test  the  model  over  a  broader  domain  in  the  phase 
space.  The  resulting  parametric  diagrams  (see  Figure  5) 
demonstrate  that  our  model  largely  reflects  the  empirical 
findings.  Boundaries  and  frequencies  are  in  good  agreement  with 
the  experimental  data,  though  in  a  few  cases,  typically  at  high 
temperatures,  the  oscillatory  domains  are  narrower  than  pre¬ 
dicted  by  the  model.  This  discrepancy  may  result,  at  least  in 
part,  from  the  sharp  drop  in  amplitude  that  occurs  near  the  Hopf 
bifurcation,  which  leads  to  oscillations  that  we  cannot  resolve 
experimentally.  Simplifications  in  the  model  that  reduce  the 
complexity  of  the  system  may  also  lead  to  quantitative 
mismatches  between  experiment  and  simulation.  The  ap¬ 
proximate  calculation  of  hydronium  ion  concentrations,  espe¬ 
cially  at  high  temperatures,  is  another  factor  likely  to  reduce 
the  accuracy  of  the  model.  Looking  at  Figure  5.  however, 
we  conclude  that  the  four-variable  model  renders  all  dynami¬ 
cal  aspects  of  the  high-temperature  reaction  rather  well. 
Consequently,  we  can  make  calculations  with  some  confi¬ 
dence  up  to  100  °C  to  further  explore  the  phase  space.  The 
dashed  lines  showr  the  likely  characteristics  of  the  reaction 
at  temperatures  near  the  normal  boiling  point  of  water.  It 
seems  probable  that,  by  decreasing  the  malonie  acid  con¬ 
centration  while  keeping  the  others  constant,  one  could  reach 
~25  Hz  oscillations  near  100  °C. 

Conclusions 

Our  experimental  results  showr  that  chemical  systems  such 
as  the  Belousov— Zhabotinsky  reaction  can  support  oscillations 
with  frequencies  above  10  Hz.  We  demonstrated  that  this 
frequency  regime  is  easily  accessible  and  studying  it  requires 
no  special  equipment.  High-frequency  oscillations  can,  however, 
be  sustained  only  for  a  few  seconds  in  closed  systems,  though 
they  should  continue  indefinitely  in  a  flow  reactor. 

Numerical  simulations  show  that  our  four-variable  model 
quantitatively  describes  all  dynamical  features  of  the  reaction 
between  40  and  80  °C.  This  allowed  us  to  map  the  phase  space 
at  higher  temperatures  and  identify  possible  pathways  leading 
to  frequencies  approaching  25  Hz.  These  experiments  would 
require  a  sealed  tube  to  avoid  evaporation,  a  stir  bar  instead  of 
a  glass  rod  immersed  in  the  reaction  mixture  from  above,  and 
perhaps  a  high-performance  temperature  control  system  able 
to  compensate  for  the  rapid  heat  production. 

In  a  high-pressure  environment,  it  might  be  possible  to  exceed 
significantly  even  the  3—4  orders  of  magnitude  increase  we 
achieved  here  over  the  typical  BZ  oscillation  frequency.  In 
principle,  the  temperature  could  be  increased  to  the  critical  point 
of  water  (374  °C)  while  maintaining  oscillations,  though  other 
factors,  including  the  increased  dissociation  of  water  and  the 
decarboxylation  of  malonie  acid,  would  have  to  be  taken  into 
account.  It  would  perhaps  be  worthwhile  to  make  calculations 
using  our  model  to  estimate  a  feasible  frequency  limit  that  could 
be  reached  experimentally  if  one  were  interested  in  increasing 
the  speed  of  chemical  computation  devices32  that  incorporate 
BZ  oscillators. 
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ABSTRACT  Many  phenomena  of  biological,  physical,  and  chemical  importance 
involve  synchronization  of  oscillatory  elements.  We  explore  here,  in  several 
geometries,  the  behavior  of  diffusively  coupled,  nanoliter  volume,  aqueous  drops 
separated  by  oil  gaps  and  containing  the  reactants  of  the  oscillatory  Belousov— 
Zhabotinsky  (BZ)  reaction.  Avariety  of  synchronous  regimes  are  found,  including 
in-  and  antiphase  oscillations,  stationary  Turing  patterns,  and  more  complex 
combinations  of  stationary  and  oscillatory  BZ  drops,  including  three-phase  pat¬ 
terns.  A  model  consisting  of  ordinary  differential  equations  based  on  a  simplified 
description  of  the  BZ  chemistry  and  diffusion  of  messenger  (primarily  inhibitory) 
species  qualitatively  reproduces  most  of  the  experimental  results. 

SECTION  Kinetics,  Spectroscopy 


time 


Studies  of  synchronized  oscillators  have  generated  sig¬ 
nificant  insights  into  a  diverse  array  of  physical,  che¬ 
mical,  and  biological  phenomena  since  the  first  known 
efforts  by  Huygens  in  the  seventeenth  century.1  In  biology, 
coupled  oscillator  approaches  have  been  applied,  for  example, 
in  modeling  the  growth  of  slime  molds,2  bimanual  coordina¬ 
tion  of  vertebrates,3-5  and  in  much  more  complex  simula¬ 
tions  of  the  human  brain  6  The  coupling  can  be  local  or  global 
(e.g.,  all-to-all)  or  a  combination,7,8  as  well  as  attractive  (ex¬ 
citatory)  or  repulsive  (inhibitory).  Global  inhibitory  coupling  in 
chemical  or  electrochemical  oscillatory  systems  led  to  the 
discovery  of  such  new  patterns  as  oscillatory  or  chaotic 
clusters9,10  Global  excitatory  coupling  usually  results  in  sys¬ 
tem-wide  in-phase  synchronization  (above  a  critical  value  of 
the  coupling  strength).11  When  the  coupling  is  purely  excita¬ 
tory  or  when  excitatory  coupling  is  dominant,12,13  a  system 
typically  has  two  states:  asynchronous,  in  which  the  oscillators 
have  random  phases,  and  synchronous,  when  all  oscillators 
are  in-phase  with  one  another.  Time  delays  associated  with 
the  coupling  can  complicate  this  simplified  picture. 

Local  inhibitory  coupling  has  been  much  less  investigated 
experimentally,  though  it  is  found  in  many  important  cases. 
Examples  of  inhibitory  coupling  include  interneuron  commu¬ 
nication  in  the  brain6  and  coupled  semiconductor  lasers.14 
Several  theoretical  works  suggest  that  inhibitory  coupling 
should  produce  many  different  synchronous  regimes  and 
multistability  between  them.15-19  Recently,  we  developed 
an  experimental  system  consisting  of  a  linear  array  of  small 
(50-300  //m  in  diameter)  identical  water  droplets  separated 
by  oil  gaps.20  Each  drop  contains  the  reactants  of  the  oscilla¬ 
tory  Belousov— Zhabotinsky  (BZ)  reaction:21,22  malonic  acid 
(MA),  bromate,  sulfuric  acid,  ferroin  (catalyst),  and  a  small 
amount  of  Ru(bpy)3,  which  serves  both  as  a  cocatalyst  and  to 
make  the  BZ  reaction  photosensitive.  Droplets  are  diffusi¬ 
vely  coupled  through  the  inhibitor,  Br2,  which  preferentially 
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dissolves  in  the  oil  gaps.  We  observed  antiphase  synchroniza¬ 
tion  and  stationary  Turing  patterns  in  this  system.20 

In  this  letter  we  report  new  synchronized  patterns  ob¬ 
served  in  two-dimensional  (2D)  arrays  of  oscillatory  BZ  micro¬ 
drops,  as  well  as  novel,  more  complex  patterns  found  in  linear 
( 1 D)  configurations  and  in  arrays  of  partially  stacked  (“  1 ,5D”) 
drops.  Thus  we  experimentally  confirm  that  inhibitory  cou¬ 
pling  is  able  to  produce  a  rich  variety  of  patterns. 

In  the  1 D  configuration,  the  most  stable  behavior  is  the 
previously  reported  antiphase  oscillation,20  in  which  each 
drop  is  1 80°  out  of  phase  with  its  neighbors.  If  we  synchronize 
the  drops  initially  with  a  flash  of  light  and  block  communica¬ 
tion  between  drops  by  adding  a  bromine  scavenger  to  the  oil 
between  drops  (see  Supporting  Information,  Methods),  then 
we  initially  observe  in-phase  synchronization,  but,  because 
there  are  small  differences  in  frequency  between  the  drops, 
this  regime  demonstrates  a  slow  phase  drift  that  soon 
destroys  the  synchrony.  If  no  scavenger  is  employed,  the 
drops  remain  in  phase  for  tens  of  oscillations,  indicating 
active  in-phase  synchronization  (see  Supporting  Information, 
Figure  SI).  Our  simulations  show,  however,  that  in-phase 
synchronization  with  inhibitory  coupling  has  a  much  smaller 
basin  of  attraction  than  antiphase  synchronization.  Experi¬ 
mentally,  delaying  the  phase  of  a  single  drop  in  a  synchro¬ 
nized  in-phase  1 D  arrangement  with  a  focused  laser  beam 
initiates  a  transition  to  an  antiphase  pattern  that  propagates 
out  from  the  delayed  drop.  Indeed,  a  quite  precise  in-phase 
initial  condition  among  all  drops  is  necessary  to  produce 
in-phase  synchronization. 

In  the  1 .5D  experiments,  when  the  drop  diameter  lies 
between  0.5d  and  d,  where  d  is  the  inner  diameter  of  the  glass 
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Figure  1.  Snapshots  (a.c.d)  and  space-time  plots  (b),  (e.  ID)  of  patterns  emerging  from  initially  in-phase  arrays  of  partially  stacked  ( 1.5D) 
BZ  drops  in  a  150//m  ID  capillary.  In  frames  b  and  e,  spikes  of  oxidation  of  ferroin  are  seen  as  short  horizontal  light  lines  (flashes)  across 
BZ-drops.  In  frames  a  and  b.  horizontal  length  of  the  frame  is  3.8  mm;  total  time  shown  in  b  is  3900  s.  In  frames  c,  d.  and  e.  horizontal  lengths 
are  1 .6.  2.2.  and  4.6  mm.  respectively.  Total  time  shown  in  e  is  1 2  000  s.  Initial  concentrations  of  reactants;  [H2S04]  =  0.08  M,  [NaBr03]  = 
0.288  M.  [MAJ  =  0.64  M  (a.  b),  0.040  M  (c-e).  [NaBrJ  =  10  mM,  [ferroin]  =  3  mM,  and  [Ru(bpy)3]  =  0.4  mM  (bpy  =  bipyridine). 


capillary,  the  drops  self-organ ize  to  form  a  zigzag  structure 
(Figure  la).  This  geometry,  in  which  each  drop  has  four 
neighbors,  two  near  and  two  more  distant,  produces  a  wealth 
of  new  stable  patterns,  some  of  which  are  shown  in  Figures  1 
and  2.  The  space— time  plots  consist  of  the  intensity  taken 
along  a  horizontal  line  through  the  center  of  the  capillary  as  a 
function  of  time.  In  Figure  lb  we  follow  the  evolution  of  60 
drops  in  one  very  complex  regime,  in  which  all  the  drops 
oscillate.  Initially,  all  drops  are  in-phase,  but  in  the  second 
cycle  shown,  drop  38  (numbering  starts  from  the  left)  experi¬ 
ences  a  delay  in  flashing.  This  phase  shift  cascades  to  the 
neighboring  drops,  generating  a  switching  wave  shown  by  the 
green  dotted  lines.  The  velocity  of  this  wave  is  about  0.8  //m/s, 
approximately  equal  to  the  linear  size  of  a  drop  divided  by  the 
period  of  oscillation.  Note,  however,  that  the  oscillation  period 
can  vary  significantly  depending  on  the  synchronous  regime. 
For  example,  the  period  of  the  initial  in-phase  regime  is  much 
shorter  than  that  of  the  final  staggered  pattern.  The  same 
behavior  emerges  in  the  simulations. 

Looking  carefully  at  Figure  lb.  we  see  that  the  drops  near 
drop  38  oscillate  180°  out-of-phase  with  their  neighbors.  A 
typical  stretch  of  out-of-phase  behavior  is  marked  by  the  blue 
oval.  Other  sets  of  neighboring  drops  assume  phase  differ¬ 
ences  close  to  ±120°,  as  highlighted  by  red  dashed  lines  in 
Figure  lb.  There  are  also  several  pairs  of  neighboring  drops 
that  oscillate  in-phase  from  time  to  time.  A  simpler  behavior  in 
1.5D,  stationary  Turing  patterns  found  at  a  lower  MA  con¬ 
centration.  is  illustrated  in  Figure  Ic.d.e.  In  the  experiment 
shown  in  Figure  1  c.d.  all  drops  initially  oscillate  in-phase  after 
a  brief  synchronizing  illumination  (Figure  lc).  Soon,  each 
drop  ceases  to  oscillate  and  reaches  a  stationary  reduced  or 
oxidized  state,  forming  the  Turing  structure  (Figure  Id). 
The  characteristic  length  of  this  pattern  is  several  drops, 
suggesting  that  the  Hiring  wavelength  can  be  larger,  perhaps 
much  larger,  than  the  distance  between  neighboring  drops. 
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The  dynamics  of  an  array  of  drops  in  a  ID  chain  that 
ultimately  reaches  a  stationary  and  more  regular  Hiring 
pattern  obtained  in  another  experiment  is  shown  in 
Figure  le.  Here,  the  characteristic  length  of  the  pattern  com¬ 
prises  exactly  two  drops  in  a  reduced  state  and  one  drop  in  an 
oxidized  state.  Note  that  before  the  system  reaches  the 
stationary  Turing  regime,  the  drops  oscillate  antiphase. 

Figure  2  demonstrates  a  novel  pattern  found  in  the  1 .5D 
geometry,  where  all  drops  in  one  row  (see  Figure  2a,d)  are 
stationary,  while  drops  in  the  other  row  oscillate  antiphase. 
Figure  2b  shows  the  initial  stage  of  these  oscillations,  and 
Figure  2c  shows  the  final  stable  regime.  These  drops  were  not 
initially  synchronized,  as  no  Ru(bpy)3  was  added  to  this 
system,  so  the  drops  begin  to  oscillate  with  random  phases 
(Figure  2b).  A  schematic  representation  of  the  pattern  beha¬ 
vior  is  shown  in  Figure  2d:  drops  in  the  bottom  row 
(represented  as  blue  or  red)  oscillate  antiphase,  while  drops 
in  the  upper  row  remain  in  a  stationary  reduced  state  (gray). 
The  spatial  period  of  this  pattern  consists  of  four  drops.  This 
pattern  is  neither  a  Turing  pattern  nor  a  simple  antiphase 
oscillatory  pattern.  It  is  difficult  to  identify  an  analogous 
behavior  in  the  homogeneous  spatially  extended  reac¬ 
tion-diffusion  system,  although  if  we  equate  stationary  drops 
to  nodes,  then  this  pattern  might  be  thought  of  as  a  discrete 
analog  of  a  standing  wave.  Another  possible  analogy  is  to  the 
oscillatory  out-of-phase  Turing  patterns  found  recently.23,24 

In  2D  geometry,  the  drops,  when  packed  at  high  density, 
spontaneously  form  a  hexagonal  array.  The  regime  consisting 
of  antiphase  oscillation  between  all  neighboring  pairs,  which 
is  extremely  stable  in  the  1  D  geometry,  is  impossible  in  a 
hexagonal  2D  geometry.  A  similar,  symmetry-generated  si¬ 
tuation  arises  in  an  analogous  Hamiltonian  system:  the  2D 
antiferromagnetic  XY  model25,26  To  resolve  this  geometrical 
constraint,  our  drops  respond  in  different  ways.  Figure  3a 
shows  a  pattern,  denoted  “jt-S”,  in  which  stationary  drops 
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Figure  2.  Pattern  consisting  of  stationary  and  antiphase  oscilla¬ 
tory  BZ  drops  in  1 .5D  geometry  (a)  Snapshot  of  capillary  with  BZ 
drops;  (b.c)  space— time  plots  (time-axis  is  vertical)  for  (b)  early 
times  during  transient  period  and  (c)  later,  after  the  stationary 
regime  is  established,  (d)  schematic  representation  of  the  pattern. 
Initial  concentrations  of  BZ  reagents  are  as  in  Figure  1  with  [MAJ  = 
0.64  M,  except  [Ru(bpy)3J  =  0.  Drop  size  is  65  /*m.  Oscillation 
period  is  177  s.  Length  of  frames  is  1.044  mm.  Times  shown  in 
frames  b  and  c  are  1068  and  798  s,  respectively. 

(black)  surrounded  by  six  oscillating  antiphase  drops  form  a 
larger  hexagonal  pattern.  Figure  3b  shows  this  pattern  sche¬ 
matically.  with  stationary  drops  in  gray,  and  blue  and  red 
drops  oscillating  out  of  phase  by  nearly  1 80°.  As  a  result,  drop 
10  (or  drops  2, 4.  7.  13.  16.  1  8.  etc.)  in  Figure  3b  is  exposed  to 
inhibitory  bromine  signals  twice  as  often  as  any  of  the  red- 
blue  pairs.  This  inhibition  is  apparently  sufficient  to  prevent  it 
from  oscillating,  pushing  it  into  a  stationary  reduced  state. 
This  2D  pattern  bears  similarities  to  the  mixed  stationary- 
oscillatory  1.5D  pattern  shown  in  Figure  2;  however,  the 
spatial  period  of  the  pattern  in  Figure  3  taken  along  any  of 
the  three  hexagonal  axes  consists  of  three  drops.  Figure  3c,d 
shows  space— time  plots  for  selected  drops,  confirming  the 
antiphase  behavior  of  neighboring  drops. 

A  second  configuration  for  drops  in  2D.  denoted  “2jt/3.” 
is  shown  in  Figure  4.  where  the  neighboring  drops  oscillate 
with  a  phase  shift  close  to  120°.  corresponding  to  the  ground 
state  of  the  2D  antiferromagnetic  XY  model.25,26  A  snapshot 
of  the  system  is  shown  in  Figure  4e,  with  the  outlined  region 
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Figure  3.  Stationary  and  antiphase  oscillatory  BZ  drops  in  2D. 
(a)  Hexagonally  packed  drops  forming  a  larger  hexagon  with  each 
stationary  drop  (black)  surrounded  by  six  oscillatory  drops.  Im¬ 
aged  area  is  0.7  x  0.7  mm.  (b)  Schematic  representation  of  the 
pattern  (stationary  drops  in  gray),  (c.d)  Space— time  plots 
(duration  2800  s)  of  drops  6.  11.  15,  14.9.  and  5  (c)  and  drops  1. 
5.  10,  15.  and  19  (d).  Concentrations  of  BZ  reagents  are  as  in 
Figure  1. 

of  interest  magnified  in  Figure  4f.  Frames  a— dare  snapshots 
of  the  region  centered  around  drop  1  in  Figure  4f  taken  after 
the  stationary  behavior  has  been  reached.  The  interval 
between  snapshots  is  about  73  s.  i.e..  1  /3  of  the  single  drop 
oscillation  period  of  220  s.  Drops  2.  4,  and  6  oxidize  in-phase 
in  frame  a.  followed  by  in-phase  oxidation  of  drop  1  and  its 
symmetrically  equivalent  counterparts  (b).  and  finally  drops 
3.  5.  and  7  (c).  The  cycle  starts  again  in  frame  (d).  Figure  4g 
shows  a  staggered  space— time  plot  of  drops  1—7,  marked 
in  (f).  This  regime  resembles  the  portions  of  Figure  1  b  with 
phase  shift  =s  120°  between  neighboring  drops.  The  most 
important  difference  between  the  conditions  that  result  in 
the  two  patterns  illustrated  in  Figure  3  (jt-S)  and  Figure  4 
(2jt/3)  is  that  jt-S  patterns  are  obtained  at  stronger  coupling, 
which  we  can  characterize  by  the  dimensionless  ratio 
DI(kL2).  Here  D  is  the  effective  diffusion  coefficient  of  the 
inhibitor  (bromine),  k  is  the  effective  (first-order)  rate  con¬ 
stant  for  the  consumption  of  the  inhibitor,  and  L  is  the 
distance  between  drops.  There  are  at  least  two  ways  to 
increase  the  coupling  strength:  (/)  by  decreasing  the  size  of 
the  water  droplets  or  (ii)  by  decreasing  the  length  of  the  oil 
gap  between  droplets.  We  explored  both  methods  (to  de¬ 
crease  the  oil  gap,  we  decreased  the  amount  of  oil  in  the 
emulsion).  With  either  method,  an  increase  in  coupling 
strength  (e.g..  by  decreasing  the  droplet  diameter  from 
1 20  to  90  /<m)  leads  to  a  transition  from  the  2jt/3  to  the 
jt-S  pattern.  More  generally,  stronger  coupling  results  in  a 
higher  proportion  of  stationary  droplets,  a  result  found  in 
both  the  experiments  and  the  simulations.  The  simulations 
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Figure  4.  Three-phase  oscillatory  clusters  in  2D.  Drop  diameter  is 
128  /<m.  Initial  concentrations:  [FLS04]  =  0.08  M,  [NaBr03]  = 
0.29  M,  [MA]  =  0.64  M,  [Ferroin]  =  3  mM.  [NaBr]  =  0.01  M,  and 
[Ru(bpy)3]  =  1 .2  mM.  Snapshots  a-d  are  taken  at  t  =  0.  773,  2  773. 
and  T,  where  7*  ==  231  s.  Snapshot  e  and  its  enlargement  (f)  show  a 
2D  configuration  of  BZ  drops,  (g)  Space— time  plot  for  drops  1  —7 
shown  in  snapshot  f;  dimensions  are  896//m  horizontal  and  2877  s 
vertical.  Period  of  oscillations  for  each  drop  in  the  space- time  plot 
is  the  same  as  T  in  snapshots  a— d. 

further  suggest  that,  in  the  limit  of  very  strong  coupling,  a 
limit  we  have  achieved  in  the  1 D  but  not  in  the  2D  experi¬ 
ments.  all  droplets  become  stationary,  resulting  in  Turing 
patterns.  At  intermediate  coupling  strengths,  the  slightly 
different  initial  conditions  may  also  contribute  to  the  selec¬ 
tion  of  different  final  patterns  by  our  system.  In  both  cases, 
we  have  patterns  consisting  of  three  groups  of  droplets 
having  either  different  phases  or  different  dynamical  beha¬ 
vior  (stationary  or  oscillatory  drops). 

To  better  understand  our  experimental  results,  we  carried 
out  computer  simulations  on  a  simplified  model  of  our  system 
SI— SI 2.  described  in  the  Supporting  Information,  Modeling. 
Because,  even  at  this  level  of  simplification,  modeling  large 
arrays  of  drops  is  computationally  taxing,  we  explored  mainly 
configurations  consisting  of  small  numbers  of  drops  to  see 
whether  the  key  features  of  our  1 D.  1 .5D,  and  2D  systems 
would  emerge. 

Typical  results  of  simulations  are  shown  in  Figure  5.  The 
regime  seen  in  Figure  2  with  stationary  and  antiphase  oscillatory 
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Figure  5.  Simulations.  (a,d)  Eight  identical  coupled  BZ  drops  for 
1 .5  D  configuration,  (a)  Gray  drops  are  stationary.  Drops  1  and  5 
oscillate  in-phase,  drops  3  and  7  also  oscillate  in-phase  but 
antiphase  to  drops  1  and  5.  (d)  Oscillations  of  drops  1  and  3, 
T  =  165.5  s,  z  is  the  oxidized  form  of  the  catalyst.  Chemical 
parameters:  h  =  [H ']  =  0.2  M,  a  =  [Br03  1  =  0.3  M,  m  =  [MA]  = 
0.5  M,  c0  (total  concentration  of  the  catalyst)  =  0.003  M.  Coupling 
parameters:  fcf  (coupling  strength  between  close  droplets  like  1 
and  2)  =  Is-1, 1/fet,  (associated  with  time  delay  in  coupling)  =  PB/ 
(partition  coefficient  for  Br2  between  water  and  oil)  =  2.5, 
pv  (volume  ratio  of  water  drops  and  oil  gaps)  =  10.  fefd  (coupling 
strength  between  distant  droplets  like  1  and  .3  in  frame  a)  =  fy/ 10, 
pv 2  =  2,fclxl  =  Pv z^fci^B.  ko  =  kfr  =  0.  (b,e)  The  minimum  number 
(three)  of  drops  simulating  the  pattern  in  Figure  3.  The  gray  drop  is 
stationary.  Drops  1  and  3  oscillate  antiphase,  as  shown  in  frame  e. 
Parameters:  h  =  0. 16  M,  a  =  0.3  M.  m  =  0.65  M.  k{  =  Is _l.  PB  = 
2.5,  pv2  =  10.  (c.f)  Twenty-eight  identical  hexagonally  packed  BZ 
drops,  (c)  Gray  drops  (7,  21.  12,  25.  16.  and  17)  are  in  the  oxidized 
SS  (oscillates  with  very  small  amplitude).  Reddish  drops  oscillate 
almost  in-phase.  Bluish  drops  also  oscillate  almost  in-phase  and 
antiphase  to  reddish  drops,  (e)  phase  shift  between  three  types 
(with  slightly  different  environments)  of  bluish  drops  (1,5,  and 
26).  Parameters:  as  in  (b.e)  except  k{  =  0.4  s  \  Other  constants 
used  in  the  FKN  model  SI— SI 2  are  specified  in  the  caption  for 
Figure  3S  in  the  Supporting  Information. 

BZ  drops  is  reproduced  in  Figure  5a.d  for  an  eight-drop  array 
in  the  1 .5D  configuration.  Stationary  Turing  patterns  re¬ 
sembling  the  experimental  patterns  shown  in  Figure  ld,e 
are  presented  in  the  Supporting  Information,  Results  (Figure 
S3).  Such  Turing  patterns  are  obtained  with  stronger  inhibi¬ 
tory  coupling  (through  bromine)  and/or  at  smaller  concen¬ 
trations  of  M  A,  as  in  the  experiments.  Note  that  a  decrease  in 
[MA]  leads  to  higher  [Br2]  in  the  drops  and  the  intervening  oil 
and  consequently  gives  stronger  inhibitory  coupling. 

The  2  D  pattern  shown  in  Figure  3  can  be  simulated  with  an 
array  of  three  linearly  coupled  drops  in  1 D  shown  in  Figure  5b,e 
(where  the  central  drop  is  nearly  stationary  and  the  two  end 
drops  oscillate  antiphase,  over  a  broad  range  of  concen¬ 
trations)  and  several  configurations  of  hexagonally  packed 
drops  shown  in  Figure  5c.  For  example,  if  we  take  only  drops 
1-7  as  labeled  in  Figure  5c,  drops  1-12.  1-16.  1-20,  or 
1—24,  we  obtain  patterns  as  in  Figure  3.  However,  if  the 
number  of  drops  in  our  cluster  is  small  (less  than  about  20), 
the  final  pattern  is  very  sensitive  to  the  addition  or  removal 
of  a  single  droplet.  Acluster  of  28  drops  (shown  in  Figure  5c) 
is  much  less  sensitive  to  this  procedure.  A  number  of 
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questions  about  the  sensitivity  of  the  pattern  to  the  number 
of  droplets,  the  coupling  strength,  the  symmetry  of  the 
droplet  configuration,  and  the  initial  conditions  remain.  Note 
that  droplets  that  oscillate  “in-phase,”  like  the  reddish  or 
bluish  droplets  in  Figure  5c,  have  a  slight  phase  shift,  as 
shown  in  Figure  5f.  An  analogous  phase  shift  can  be  seen  in 
all  our  experimental  results.  The  origin  of  this  shift  in  the 
simulations  is  the  different  number  of  neighbors  for  bound¬ 
ary  and  “inner”  droplets. 

The  three-phase  oscillatory  patterns  shown  in  Figure  4  can 
be  simulated  with  just  three  circularly  coupled  BZ  drops  plus 
three  oil  drops  (like  drops  1,2,  and  7  in  Figure  5c,  for  example; 
see  also  Supporting  Information,  Results,  Figure  S4).  Curi¬ 
ously,  the  three-phase  pattern  of  Figure  3  cannot  be  modeled 
with  three  circularly  coupled  drops.  Both  of  the  two  three- 
phase  patterns  in  2D,  jt-S  and  2jt/3,  are  found  in  a  cluster  of  six 
identical  drops  connected  in  3D  to  form  an  octahedron 
(Supporting  Information,  Results,  Figure  S5).  This  unique 
configuration  allows  us  to  have  periodic  boundary  conditions, 
an  even  number  of  neighbors  (4  vs  6  in  hexagons),  and  a  total 
number  of  drops  divisible  by  3  (the  spatial  period  for  these 
two  patterns).  Simulations  reveal  that  an  increase  in  coupling 
strength  leads  to  a  transition  from  2jt/3  to  jt-S  patterns,  just  as 
we  observe  in  our  experiments.  However,  there  is  no  clear 
physical  correspondence  between  the  experimental  configu¬ 
ration  of  hexagonally  packed  drops  and  the  simulations  on 
3D  octahedrons. 

In  our  simulations,  we  used  the  well-known  FKN  mecha¬ 
nism  for  the  BZ  reaction,27  in  which  Br02*  radical  serves  as  an 
excitatory  messenger  between  droplets,  and  Br2  acts  as  an 
inhibitory  messenger.  Excluding  Br02*  as  a  messenger  in  the 
model  causes  almost  no  change  in  the  calculated  patterns, 
which  suggests  that  this  species  plays  practically  no  role  in  the 
coupling. 

The  formation  of  two-phase  and  three-phase  structures, 
i.e.,  sets  of  drops  that  oscillate  with  a  common  phase  (or  are 
stationary)  in  our  system,  may  be  viewed  as  a  generalization 
of  the  chemical  quorum  sensing  recently  reported1 3  in  sets  of 
particles  loaded  with  catalyst  and  immersed  in  a  BZ  solution. 
Here,  too,  the  individual  drops  behave  independently  until 
brought  together  so  that  their  density  (or  coupling  strength) 
exceeds  a  critical  value.  They  then  begin  to  exhibit  more 
complex  dynamical  behavior,  in  our  case  antiphase  and  three- 
phase  patterns  as  well  as  the  in-phase  oscillation  observed  in 
the  particle  experiments.  Our  oscillatory  structures  also  re¬ 
semble  the  clusters  observed  in  experiments  on  globally 
coupled  electrochemical  oscillators.10  The  larger  the  coupling 
strength,  the  larger  the  number  of  different  synchronous 
groups  (up  to  the  limit  of  Turing  patterns).  Therefore  in 
addition  to  the  “quorum  sensing”,  our  system  exhibits  a 
primitive  kind  of  “differentiation”. 

Since  the  coupling  of  our  oscillatory  BZ  drops  is  primarily 
inhibitory,  and  they  can  form  stationary  Turing-like  patterns  (if 
the  coupling  is  strong  enough),  this  system  can  be  viewed  as  a 
discrete  manifestation  of  the  Turing— Hopf  interaction  that 
generates  many  interesting  patterns  in  continuous  reaction- 
diffusion  systems,  such  as  oscillatory  in-phase  and  antiphase 
Turing  patterns  (the  latter  is  analogous  to  our  antiphase 
regime).23,24 
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All  configurations  that  we  have  studied  computationally 
demonstrate  multiple  synchronous  regimes,  i.e.,  regimes  in 
which  each  drop  has  the  same  total  period  (though  possibly 
with  different  numbers  of  maxima  during  this  period).  We 
expect  to  find  some  of  these  regimes  in  further  experiments 
on  our  system  of  coupled  BZ  droplets,  which  mimics  diffusive 
inhibitory  coupling  among  biological  objects,  a  situation  that 
is  ubiquitous  in  nature. 

SUPPORTING  INFORMATION  AVAILABLE  Methods,  mode¬ 
ling,  and  results  as  Figures  SI  — S5.  This  material  is  available  free  of 
charge  via  the  Internet  at  http://pubs.acs.org. 
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